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AND MIXED MODE FRACTURE: A SUBLAMINAI E/PLY LEVEL ANALYSIS 

AND A COMPUTER CODE 

R.R. Vallsetty* and C.C. Chamls 
National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


SUMMARY 

A computer code Is presented for the sublamlnate/ply level analysis of 
composite structures. This code Is useful for obtaining stresses In regions 
affected by delamlnatlons , transverse cracks, and discontinuities related to 
Inherent fabrication anomalies, geometric configurations, and loading condi- 
tions. Particular attention Is focussed on those layers or groups of layers 
( sublamlnates) which are Immediately affected by the Inherent flaws. These 
layers are analyzed as homogeneous bodies In equilibrium and In Isolation from 
the rest of the laminate. The theoretical model used to analyze the Individual 
layers allows the relevant stresses and displacements near discontinuities to 
be represented In the form of pure exponential -decay- type functions which are 
selected to eliminate the exponential-precision-related difficulties In sub- 
lamlnate/ply level analysis. Thus, sublamlnate analysis can be conducted with- 
out any restriction on the maximum number of layers, delamlnatlons, transverse 
cracks, or other types of discontinuities. In conjunction with the strain 
energy release rate (SERR) concept and composite micromechanics, this computa 
tlonal procedure Is used to model select cases of end-notch and mixed mode 
fracture specimens. The computed stresses are In good agreement with those 
from a three-dimensional finite element analysis. Also SERRs compare well with 
limited available experimental data. 


INTRODUCTION 

Laminated composite structures exhibit a number of different failure 
modes. These Include fiber-matrix debonding within Individual layers, delaml 
nation or separation of adjacent layers, transverse cracking through one or 
more layers, and fiber fracture. Often Interlaminar delamlnatlon Initiates at 
sites where complex states of stress with steep gradients exist. These sites 
Include regions near free edges, ply terminations, cutouts, voids, holes Inad- 
vertantly damaged areas, and defects resulting from fabrication processes. As 
a result, an accurate assessment of the Interlaminar fracture toughness param 
eters and stress states which initiate delamlnatlon Is of critical Importance 
In a composite design. 

Several methods are available to determine fracture toughness and can be 
classified into the following groups: experimental, numerical, and analytical, 

lo begin with, the experimental procedures used to determine fracture toughness 
are edge delamination, double- cantl lever- beam, cracked- lap- shear , biaxial 
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Inter laminar fracture, and three point bend tests for both Mode I (end notch 
fracture, tNt ) and mixed mode fracture (MMF) Including Mode 1 and Mode II. 
because of the test simplicity and the ability to delineate between Mode II and 
mixed mode fracture, the three-point bend tests are now being Investigated at 
NASA Lewis Research Center. Results for Interlaminar fracture toughness of 
unidirectional graphite/epoxy composites are presented In reference 1. 

Next, there are numerical solutions. In general, these are developed 
using either finite difference, finite element, finite-element boundary 
Integral concepts, boundary layer, energy based approximations, and extended 
Galerkln procedures. They require extensive Idealizations to model the 
response variables In the critical regions and the numerical accuracy of the 
results Is not always assured. Issues related to this topic are discussed In 
references 2 to 5. 


Conventional laminate theories form the core of the final group, the ana- 
lytical approaches. These are best suited for obtaining the total strain 
energy release rates. Fracture related Interlaminar stresses which are neces 
sary for delineating between different modes of crack propagation cannot be 
obtained from these theories. This Is because the theories are based on plane 
stress and global displacement assumptions. The assumption of plane stress 
precludes the possibility of considering the Interlaminar stresses at the out 
set of the formulation. 


This led to the development of ply/sublamlnate level analysis (refs. 2 
and 4) of composite laminates. Here, layers or selected groups of layers 
(treating each group as a homogeneous body) are analyzed with the assumption 
that they are In equilibrium and Independent of each other. Thus, the Inter- 
laminar stresses are Introduced In the formulation at an earlier stage. 
Initially, these stresses are assumed to be unknown. Then considering the 
equilibrium of each layer/sublamlnate and compatibility at each Interlaminar 
surface, a set of coupled equations In terms of global deformation variables 
and Interlaminar stresses Is generated. Subsequent to a successful solution, 
the fracture toughness and SERRs are easily computed. This approach was Ini- 
tially applied to the estimation of free-edge stresses (ref. 2), and later, 
successfully extended to both an edge delamlnatlon problem (ref. 6) and the 
analysis of compression and double-lap fracture specimens (ref. 7). 


This report describes the Implementation of sublamlnate analysis to 
general composite laminates. This analytical technique has been Incorporated 
Into a computer code which Is listed In appendix I. Ihe code Is written with 
flexibility In regards to the types of materials and loads, number of layers, 
and boundary conditions as well as cracks which may be In the plane of the 
laminate or through the thickness. The versatility of this method and the 
numerical economy In conducting design related parametric studies are demon 
strated through the simulation of three-point bend tests of ENF and MMF speci- 
mens. A brief description of sublamlnate analysis and ply level equations are 
presented and an Illustration of the basic steps Involved In ply-modeling a 
given structure and setting up the governing equations Is Included. 


SUBLAM1NA1 E/PLY LEVEL ANALYSIS 


This 
layers In 


analysis which operates on a ply by ply basis may be applied to all 
a laminate or only to those few layers which surround defects and 
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where the stresses are adversely affected. The remaining layers can be grouped 
Into sublamlnates and may be treated as homogeneous layers. Interlaminar 
stresses are Initially assumed to be unknown. Enforcement of displacement and 
traction continuity at laminar Interfaces leads to the final equations, which 
are subsequently solved. Thus, a better engineering approximation Is achieved 
concerning the interlaminar stresses. The success of the approach, then, 
depends on the engineering model that will be used to analyze each layer. The 
model should be capable of securing " layer equilibrium " as discussed below. 

Just as It is possible to combine different layers Into a single unit or 
a sublaminate and treat the resulting unit as a single layer. It Is also pos- 
sible to Idealize a single layer as composed of a finite number of thin sub- 
layers. For the purpose of sublaminate/ply level analysis, these sublayers 
once again can be treated as Individual units. This situation Is, therefore, 
similar to employing either a coarse mesh or a fine mesh at selective locations 
through the thickness In a finite element analysis. The ensuing results for 
Interlaminar stresses depend upon the Idealization In terms of sublamlnates, 
layers, and sublayers and to what extent they are affected Is discussed In 
reference 4. 

To achieve "layer equilibrium," the model must be capable of dealing with 
prescribed stress/displacement conditions on transverse Interlaminar surfaces 
and afford specification of at least five boundary conditions (three force 
resultants or deflections and two moment resultants or rotations) per layer 
edge. Two such models are available now. The first one, developed by Pagano 
(ref. 2), is a higher order theory based on Relssner's variational principle. 
Ihe second one, developed by Vallsetty (ref. 8) from an Iterative procedure, 
is a refined homogeneous plate bending theory. The later model was selected 
to analyze the layers in the present program because It provides explicit dis- 
tributions through the thickness for stresses and displacements. 

A summary of the basic equations for the generic ply shown In figure 1 Is 
given in appendix II. The equations are deduced from the homogeneous plate 
theory (ref. 8) for the case of cylindrical bending. The basis for these equa- 
tions and their derivation are presented In references 6 and 8. 

Suppose the analysis Is to be applied to N layers or a combination of 
layers and sublamlnates. Ihe overall equations of equilibrium and the constl 
tutlve relationships will constitute a set of 8N equations In terms of like 
variables (2N displacements, N rotations, 3N force resultants, and 2N 
moment resultants). This set Is further supplemented by additional 2(N 1) 
equations, which are necessary for the simultaneous solution of 2( N- 1 ) Inter- 
laminar stresses, when the displacement continuity Is enforced at the (N 1) 
Interlaminar surfaces. The transverse shear and normal stresses at the Inter- 
laminar surfaces, If they are not prescribed explicitly (for example, as on the 
top and bottom surfaces of a laminate or on the cracked Interlaminar surfaces), 
are treated as unknowns. From among this set of equations the force and moment 
resultant variables can be eliminated with the aid of constitutive relations. 
This leaves a set of (5N 2) coupled ordinary differential equations to be 
solved for 2N displacement variables, N rotations, and 2(N-1) Interlaminar 
stresses . 

These observations are based upon two requirements: (1) the laminate is 

In generalized plane strain, and (2) there Is no cross-sectional warping. The 
latter condition Is satisfied when cross- ply laminates and symmetric laminates 
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are used. However, for angle ply laminates and laminates of more general con 
struction, the warping cannot be neglected. In this case, the number of varl 
ables Increase and likewise the number of equations to be solved. 

Ihe final system of differential equations can be represented In the 
following form: 

4 4N- 2 

£ £ SK^ } f k) = p i; 1 = 1,2 4N-2 (1) 

k*0 j=l 

where p represents the applied loading, f Is the vector of governing 
variables (2N average displacements and 2N-2 Interlaminar stresses), and 
the superscript k Is the order of the In-plane derivative. The coefficient 
matrices ( [ SK ] ) depend on layer thickness and material properties. The total 
solution Is constructed by obtaining: (1) a particular solution corresponding 

to the applied load, and (2) other solutions from the homogeneous part of the 
differential equations. The analysis shows that there are (8N-2) roots or 
eigenvalues, therefore an equal number of homogeneous or eigen solutions. An 
Inspection of the eigenvalues reveals the nature of these solutions. There are 
six zero roots and the remaining are either real or complex. The solution cor- 
responding to the zero roots is Identified to be of shear deformation type and 
the remaining are of local type which are prominent only in the edge regions 
of the laminate. 

There are (8N-2) Integration constants which require specification of a 
like number of boundary conditions at the laminate edges. These are estab- 
lished on the basis of physical reasoning. Since global laminate equilibrium 
is translated into local layer equilibrium In this analysis, boundary condi- 
tions are to be specified at the layer edges insofar as possible. Among such 
conditions, the engineering boundary conditions appear to be a natural choice. 
Consequently, one can prescribe either force resultant or an average In plane 
displacement, moment resultant or an average rotation and shear resultant or 
an average transverse displacement at each layer edge. The remaining ( 2 N - 2 ) 
conditions are specified In terms of either interlaminar transverse shear 
stress or its in plane derivative at the In-plane edges of interlaminar 
surfaces . 

By retaining the local solutions at only one laminate edge, the number of 
boundary conditions can be reduced to (4N*2). The edge at which the local 
solutions are to be retained depends on the particular problem In question. 

This procedure does not threaten the laminate equilibrium as it can be secured 
from the global perspective with the aid of a shear deformation type solution 
(which Is controlled by the six zero roots) alone. However, for this procedure 
to be successful, the layer-defined force and moment resultants at the edge, 
where the local solutions are dropped, must be specified to be statically 
equivalent to the local total loading applied to the whole laminate edge. At 
the edge where the local solutions are retained the boundary conditions can 
still be specified on each individual layer. 

This maneuver makes the sublaminate/ply level analysis powerful, for it 
allows equilibrium in any number of layers to be considered independently. 
Usually, the number of layers to which this type of analysis Is applied is 
restricted by a computer precision requirement. This requirement puts a limit 
on the magnitude of the real roots and real part of the complex roots. When 
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the local solutions are obtained for each laminate edge, one at a time, the 
mode shape of the eigen solutions can be constructed In the form of pure decay 
type solutions which are valid only at one edge. Thus the exponential- 
function related precision problem are effectively eliminated. This Is a sig- 
nificant contribution, because now sublaminate/ply level analysis can be 
extended to laminates with discontinuities In loads, materials, and geometry 
along the axis of the laminate and without any restriction as to the number of 
layers . 

For cases Involving axial discontinuities, the sublamlnate analysis pro 
cedure essentially remains the same. The first step is to divide the structure 
Into several regions or subunits such that each subunit Is free of any discon- 
tinuity related to loads, materials, or geometry. Each subunit, then, will be 
in the form of a laminate loaded by continuously varying loads. Materials may 
differ for layers within the subunits. The concentrated loads on the original 
structure translate into shear loads applied to the edges of the subunits. The 
delaminations can be accounted for by treating the separate portions of the 
laminate as different sublaminates abutting the uncracked laminate which may be 
treated as just another sublamlnate. This process Is described In detail In 
the next section with an example of ENF and MMF specimens. 

Formal solutions In terms of Integration constants can then be constructed 
for each subunit using the aforementioned procedure. The postulation of bound- 
ary conditions and continuity conditions between different subunits provides 
equations for the solution of Integration constants. Typical continuity condi- 
tions Involve specifying axial force resultant, shear resultant, moment result 
ant, axial displacement, transverse displacement, and cross-sectional rotation 
for each layer and transverse shear stress and Its derivative for each Inter 
laminar surface. 

lhe computer code containing these procedures Is presented In appendix I. 
Instructions for the Input and output are given In appendix 111. Typical Input 
and output for the ENF specimen are given In appendices IV and V, respectively. 


IDEALIZATION OF THREt POINT BEND TESTS 

The sublaminate/ply level analysis Is now applied to estimate the Inter- 
laminar fracture stresses In the three- point bend test specimens. Figure 2 
presents the schematics for the three-point bend tests of the ENF and MMF 
specimens. For the MMF specimen, the right support Is provided on the bottom 
surface of the top flange Instead of the bottom flange. The laminate is made 
of 24 graphite/epoxy (AS/E) layers and five interlaminar resin layers as shown 
in figure 3. The manufacturing process gives rise to the resin rich region 
between adjacent fiber-matrix layers whose thickness usually amounts to a tenth 
of the fiber diameter and In general depends on the fiber volume ratio ( F VR ) . 
This region is referred to as the Interply resin layer and is considered in 
order to accurately compute those stresses which induce Interply delamination. 
Two different FVRs, 0.3 and 0.6, are considered In the computations. At these 
FVRs, the Interply resin layer thicknesses are 0.0001854 In. and 0.00004323 In., 
respectively. The graphite/epoxy layers have the following properties: 
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At tVK -= 0.3 


El -- 11-96 ms 1 , E-| = 0.921 msl , E] | = 1.01 msl , G 7 \ = 0.245 msl , 

G |_ | = 0.368 msl, i>l 1 - 0.348, vj| - 0.546 
ply thickness = 0.005 In. 

At F VR = 0.6 

El = 21.15 msl, E-j = 1.238 msl, Eg j = 1.323 msl, G 77 ^ 0.362 msl, 

Gli - 0.623 msl, «L 1 - 0.282, vgy - 0.423 
ply thickness = 0.005 In. 

where E Is the elastic modulus, G Is the shear modulus and v Is Poisson's 
ratios, L represents a direction along the fibers and 1, a direction, normal 
to the fibers. 


As Inspection reveals that In order to apply sublamlnate/ply level anal- 
ysis the specimen (fig. 2 ) must be divided Into the following four regions: 


Region 1 : 0 
Region 2: L/2 
Region 3: (L-a) 
Region 4: (L-a) 


< x 2 < L/2; - H/2 < z < H/2 , 

< x 2 < ( L- 1 ) ; -H/2 < z < H/2, 

< x 2 < L; 0 < z < H/2, 

< x 2 < L; -H/2 < z < 0. 


( 2 ) 


where a, L and H denote the length of the crack, the length of the specimen 
and its thickness, respectively. The length of the specimen Is 4 In. and 
(x 2 , z) Is a global coordinate system. Region 1 Is Identical to Region 2 and 
Regions 3 and 4 are Identical. Now, the analysis can be applied to each of 
these four regions obtaining four sets of formal solutions In terms of Integra 
tlon constants. The total number of these constants will be 680. 


Before proceeding to setup such a solution. It should be asked whether the 
analysis can be simplified any further. The answer will be obvious. If It Is 
recognized that the objects of Investigation are the stresses In the Interlami- 
nar resin layer in the region ahead of the crack. To obtain these stresses, It 
is not necessary to apply the analysis to all of the layers within each of the 
four regions. Recognizing this, It can be assumed that the 12 graphite/epoxy 
and two interlaminar resin layers (fig. 3 ) can be treated as single units or 
sublaminates in the third and fourth regions and in the top and bottom halves 
of the first and second regions. These sublaminates can further be treated as 
homogeneous orthotropic plates represented by their effective moduli. The 
resulting Idealization Is shown in figure 4. 

Now, the analysis Is applied only to the first and second regions. Since 
the sublaminates In the third and fourth regions are treated as homogeneous 
plates, the governing solutions for these regions, which are of the shear 
deformation type, can be obtained directly. The order of the resulting set of 
all equations will be 56. This can further be reduced to 32 by retaining only 
the local solution near the crack In Region 2. This solution yields the frac- 
ture related stresses In the resin layer ahead of the crack. The other local 
solutions, which are suppressed, control the stress variation near the concen- 
trated load and at the left support. 
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Equations for Solution 


The solution procedure for the specimen Idealized In figure 4 Is described 
here. For this purpose, the laminate in Region 1 Is selected. Let "1," "2," 
and "3" denote, respectively, the top, middle, and bottom sublamlnates . The 
outer sublamlnates each contain 12 graphite/epoxy and two resin layers. The 
middle one contains just one Interply resin layer. For the analysis, these 
sublamlnates are treated as homogeneous plies. It Is convenient to use local 
coordinate systems, (X 2 ,z; k+1,2; and 3), for presenting stress and displace- 
ment distributions. The ply middle surfaces are represented by Z = 0 and c 
denotes ply semithicknesses . 

For each layer, stress and displacement distributions through the thick- 
ness are available. These were discussed earlier in appendix II along with 
relevant equilibrium and constitutive equations In terms of engineering kine- 
matic and force variables and Interlaminar stresses. As shown In figure 1, 
these variables are meaningful for only those plies for which they are defined. 
When It Is necessary to Indicate such ply dependency, these variables are 
superscripted with a ply Identification number. 


Flereln, let U 2 and W be the components In the X 2 and z directions, 

respectively of the displacement of the z = 0 surface. Also, let a 2 Z and 

a be the transverse shear and normal stresses. Using this notation, the 

ZZ 1231232 

unknown variables can be listed as U 2 , U 2 , U 2> W , W , W , <, 2z ( x 2 , z = c), 

o 2z ( x 2 » z = c ). °z z ( x 2 ,z = c )’ and °zz^ x 2 ,z = c )* * n *^ ese variables 

constitute a set of ten unknowns. An Identical set of variables also govern 
the deformation in Region 2. For sublamlnates defined In the third and fourth 
regions, U 2 and W are the only governing variables. 

Equations 1-2 and 1-3 of appendix II may be combined to produce 


M 2,22 + C m 2,2 + q 2 = ° 


(3) 


where 




m 2 

moment resultant associated with 

x 2 

direction. 

c 

ply semi thickness , 




m 2 

°2z( x 2» c ) + °2z( x 2* ~ c ). and 




d2 

o zz (x 2 , c) - o zz (x 2 , -c). 




The 

the 

equilibrium conditions (eq. (3), and eq. 
three plies and the following continuity 

(1-1) of appendix II) for each of 
requirements: 
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4 = u 3 ; 0 < x ? < L/2 , z 1 = -c-| , 

2 3 

u 2 = u ? ; 0 < x ? < L/2, z 2 = - c 2 , 

(4) 

1 2 

w = w ; 0 < x 2 < L/2, z^ = -c^ , and 
2 3 

w = w ; 0 < x 2 < L/2, z 2 = - c 2> 

for the displacements at the two laminar Interfaces make up the desired set of 
equations for Region 1, where U 2 . and w are the components of the dis- 
placement In the X 2 and z directions. Equation 4 can be expressed solely 
In terms of the generalized unknown variables of the problem with the aid of 
tqs. 1-16 and 1-17 of appendix II which are thlcknesswlse displacement distri- 
butions. An Identical set of equations follows for the second region. The 
cracked surfaces In Regions 3 and 4 are free of the transverse stresses. 
Therefore, the solution for the sublamlnates In these regions degenerates Into 
a shear deformation type solution with a shear correction factor of 1. 

These four sets of equations are solved separately In a straightforward 
manner and formal solutions In terms of Integration constants are constructed. 
There are 56 constants which may be determined through the boundary conditions 
at X 2 = 0 and X 2 = L and the continuity conditions between the first 
and second regions at X 2 = L/2, between the second and third regions at 
X 2 n (L-a), and between the second and fourth regions at X 2 = (L-a). For 
presenting these conditions, superscripts "4," "5," and "6" are used to denote 
the top, middle, and bottom layers In Region 2 and "7" and "8" are for the sub 
laminates in the third and fourth regions. As mentioned before, superscripts 
"1," "2," and "3" apply similarly to the variables In the first region. 

Boundary conditions at X 2 - 0: 


N 2 

= o, 

M 2 = 

o, q 2 = 

0 > <*2 £ ( *2 * 

i 

o 

ii 

O 

2 

= o, 

2 

2 

2 , 

-c 2 ) = 0 

N 2 

M 2 * 

o, q 2 = 

0 > p z ^2 * 

3 

U 2 

(x 2 , 

- C 3 } 

- 0, w 3 

( X 2 * “ c 2^ ~ 

0 M 3 = 0 
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Continuity conditions at x 2 - L/2: 


1 

4 i 

4 1 

4 

1 

4 

N 2 

- n 2 , m 2 

= m 2 , q 2 

- + c *2z^ X 2* ^ 

1 > - °2z x 2 • - C 4> 

1 

4 1 

4 . .1 

,4 

1 , 

4 , 

U 2 

a U 2’ *2 

= 4>2* w 

* w , 

a ?z X 2’ " c l ) 

~ ° 2 z x 2 9 C 4 9 2 

2 

5 2 

..8 r 2 


2 , 

5 , 

N 2 

= n 2 , m 2 

-- m 2 , Q 2 

= Q 2' 

°2z X 2 ’ ~ C 2 

- z ^2 9 _C 2 

, 2 

5 2 

8 ,,2 

,,5 

2 , 

5 , 

U 2 

- U 2 , 4>2 

= <t> 2 . w 

= w , 

c 2 

~ a 2 z( x 2 9 ~ C 2 9 2 


- n 2 , m 8 

= M 8 , Q 8 

= Q2+6O, 


.,3 

.,8 3 

8 , ,3 

. .6 



U 2 

- U 2 , <t> 2 

- 4> 2 » w 

= w , 




Continuity 

conditions 

at 

x 2 = (L-a) 

• 





fcl 4 

..7 ..4 

..7 

4 

7 



N 2 

= n 2 , m 2 = 

^ 2 . 

^2 

- Q 2 



,,4 

,,7 4 

7 

. ,4 

„ 7 



U 2 

= U 2 1 <t>2 ~ 

4>2» 

W 

= W 

Continuity 

conditions 

at 

*2 = (L-a) 

: 





..6 

..8 .,6 

,.8 


R 



N 2 

- ^ 2 * ^ 2 ~ 

^ 2 » 

Q 2 

= Q 2 



,,6 

,,8 6 

8 

, ,6 

, ,8 



U 2 

' U 2’ ^2 

V 

W 

= W 


Boundary conditions at x 2 = ( L- a ) : 

N 8 = 0, M 8 = 0, Q 8 = 0, 

4 6 

®2 Z ( x 2 ’ “ c 4^ = 0, °2z^ X 2 ’ “ c 6^ = ^ 

Boundary conditions at x 2 •= L for ENF specimen: 

= 0 , m ' = 0 , Ql -- 0 , 

N® - 0, M 8 = 0, w 8 (x 2 , -c 8 ) = 0 

Boundary conditions at X 2 - L for MMF specimen: 

n 2 = 0, m' = 0, w / (x 2 , -c 7 ) = 0, 

N 8 = 0, M 8 - 0, Q 8 - 0. 
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( 6 ) 


(7) 


( 8 ) 


(9) 


( 10 ) 


( 11 ) 


whore (k = 1,2,. ..8) are ply semi thl cknesses , dp are cross-sectional 

rotations, Np are force resultants and Mp are moment resultants. 

The applied load of 120 lb Is represented as a jump In shear resultants In 
the outer layers In equation 6. It Is assumed that the central Interlaminar 
resin layer Is too thin to carry any net shear force and the shear stress In 
the outer sublamlnates away from the point of application of the load Is uni- 
form through the thickness. Similarly, the conditions on the transverse shear 
stress In equation 6 are arrived at under the assumption that the points at 
which the stress Is prescribed are considerably away from the point of applica- 
tion of the load; therefore, the stress Is unaffected. Additionally, the crack 
Is modeled as having a finite thickness, the thickness being that of the Inter- 
ply resin layer. The Interply resin layer Is destroyed as the crack advances 
so that the crack-tip (Idealized as having finite thickness) Is stressfree. 

This Is Indicated by equation 9. 

The approach discussed so far makes use of the complete form of the eigen 
solutions to represent the response variable In the first two regions. This 
yields stress and displacement distributions which show not only the effect of 
the crack-tip but also that of the concentrated load at the left support. 

Since the objective Is the crack-tip stress field, the form of eigen solutions 
can be modified to include only the decay type solutions at the crack-tip. 

This degenerates formal representation of the solution In Region 1 into that 
of a shear deformation theory, as is also the case with the third and fourth 
regions. As a result, the total number of Integration constants and conse- 
quently the number of boundary/continuity conditions are reduced to 36. The 
conditions at xp = 0 and L/2 are to be modified as shown below: 

Boundary conditions at x? = 0: 

U 2 ( X 2 * °) = 0 * 

u^( *2 , 0) - 0 , (12) 

t hJ t nJ x Ci - X c 3 = 0 
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Continuity conditions at X 2 = L/2: 

ul (* 2 • °) = u 2< x 2 ’ °) 

w 2 (x 2 , 0) = w 5 (x 2 , 0) 

2 5 

( *2 9 = d 2^ X 2 9 ^ (13 

N 2 + N 2 + N 2 t N 2 t N 2 t N 2 

q] + «■ * <4 + Q 2 * Q 2 + 120 

M 2 f M 2 * M 2 * N 2 x c i ' N 2 x c 3 = + M 2 + x c 4 - x c & 

The procedure described so far is programmed Into a general purpose code. 
A description of this code (appendix I) together with instructions for input 
and output (appendix III) are Included. 


RESULTS AND DISCUSSION 

Figures 5 and 6 display the stress field in the interply layer ahead of 
the crack-tip for the ENF and MMF specimens, respectively. The crack length 
is 1 in. and the fiber volume ration (FVR) is 0.6. The ENF specimen yields 
momotoni cal ly decreasing compressive normal stresses, whereas the MMF specimen 
produces tensile stresses. For both specimens, the transverse shear stress is 
positive near the crack-tip. Within a distance of 1.5 H from the crack- tip, 
all stresses approach engineering values. The compressive nature of the trans- 
verse normal stress for the ENF specimen Indicates that this is a pure Mode II 
crack. For the MMF specimen, as anticipated, both transverse shear and normal 
stresses are tensile. Figures 7 and 8 show the behavior of peak stresses as 
the crack advances. Both types of specimens show monotonically increasing 
crack tip stresses. 

Strain energy release rates (SERR) are computed by monitoring changes in 
the total compliance of the structure. Figure 9, which shows the total energy 
release rate plotted against the crack length, indicates a stable crack growth. 
Both the ENF and MMF specimens predict the same total SERR. This result was 
also anticipated since, in the case of the MMF specimen, the unsupported 
delaminated flange will not pick up any load. However, in the case of the ENF 
specimen, the possibility for load transfer exists as both delaminated flanges 
come into contact as the applied load Increases. Since this contact is not 
incorporated into the analysis, the same results are predicted. 

In order to study this progressive contact, the relative displacement 
between the cracked surfaces should be restricted so that cracked sublaminates 
do not run into one another. If they do the problem becomes nonlinear and is 
beyond the scope of the present investigation. However, a simple idealization 
of the practical situation using an ENF specimen assumes complete contact 
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between the cracked surfaces and that the resulting contact stresses are negll 
glble. Then, the condition on the transverse displacement In equation 10, 
accordingly, Is replaced by the following: 

7 8 

w (x 2 , -c 7 ) = w (x 2 , Cg) (14) 

Ihls model Is named ENF- II. Results for the stress field ahead of the crack- 
tip and total SERR using this model are presented In figures 10 and 11, respec- 
tively. The FVR for these results Is 0.6. Both normal stresses are zero and 
In figure 10 the shear stress follows the previously observed trend. In this 
case, the crack plane becomes a plane of symmetry. This does not permit the 
normal stresses to develop on the crack plane. This total SERR Is only about 
one-third of that which Is available In the first two models. As shown In 
figure 11, these results compare well with engineering theory which are 
slightly conservative. Reference 1 presents some experimental data on the 
total SERR for the ENF specimen. The range of the total SERR shown In 
figure 11 matches with the experimentally measured range of 3.0 to 3.6 psl-ln. 
The role of the local stress field on crack propagation can be established by 
measuring Individual mode contributions, but that study Is deferred for a 
future Investigation. 

In order to assess the present study, the Interply stresses ahead of the 
crack tip are compared In figure 12 with the three-dimensional finite element 
results of reference 1 for the ENF specimen. With the exception of small dif- 
ferences In peak stresses at the crack-tip, both analyses provide the same 
trends. The finite element results for the transverse shear stress exhibit an 
oscillatory behavior. This was attributed. In reference 1, to the particular 
fashion In which the finite element grid ahead of the crack-tip was updated as 
the crack length Increased. 


CONCLUSIONS 

The analysis presented herein can be extended easily to cracks through the 
thickness. In this case the eigen solutions are to be computed each time the 
crack extends. Discontinuities related to geometry and loads can also be han- 
dled In an analogous manner. Similarly, flexibility Is available to consider 
delamlnatlons anywhere through the thickness and along the axis of the lami- 
nate. These features make the present analysis Ideally suited for design 
related parametric studies when: (1) different configurations and the effect 

of materials are to be evaluated, and (2) fracture characteristics need to be 
established. Ihe following conclusions are reached based on this study: 

1. After the eigen solutions are established In each region, the order of 
the final system of equations representing boundary/symmetry/continuity condi- 
tions Is reduced to 36 for the specimen Idealizations. This results In an ana- 
lysis which Is highly economical. 

2. The ability to represent critical local stress fields In terms of pure 
decay type solutions eliminates the exponential - f unction related precision 
problems. This effectively removes the cap on the number of layers that can 
be studied on a ply-by-ply basis. 
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3. Different types of restraints can be handled routinely, as Indicated 

by ENF , ENF II and MMF models, without Incurring any major changes In the final 
systems of equations. 

4. The models also Illustrate how the stress field ahead of the crack tip 
responds to simple changes In the edge supports. This knowledge Is Important 
In designing fracture specimens for different modes of crack propagation. 
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APPENDIX I: SOURCE CODE 


OOOllOO 

0001200 

0001300 

0001400 

0001500 

0001600 

0001700 

0001800 

0001900 

0002000 

0002100 

0002200 

0002300 

0002400 

0002500 

0002600 

0002700 

0002800 

0002900 

0003000 

0003100 

0003200 

0003300 

0003400 

0003500 

0003600 

0003700 

0003800 

0003900 

0004000 

0004100 

0004200 

0004300 

0004400 

0004500 

0004600 

0004700 

0004800 

0004900 

0005000 

0005100 

0005200 

0005300 

0005400 

0005500 

0005600 

0005700 

0005800 

0005900 

0006000 

0006100 

0006200 

0006300 

0006400 


PROGRAM WAIN 
C CDABS=CABS 

C DCMPLX REWOVE THIS DECLARATION 

C REAL*8 

C REAL*4 RCW 

C COMPLEX* 16 


c 

Ell 

S32 

CB16 

KQ1 

CURVTR 

2 

NZR 

12 

CRTB 

22 

c 

E22 

S33 

CB22 

KQ2 

NS REG 

3 

NRR 

13 

CRTC 

23 

c 

E33 

S44 

CB26 

KW1 

NWT 

4 

NCR 

14 

ISGN 

24 

c 

NU12 

S45 

CB66 

KM2 

NSL 

5 

NBCR 

15 



c 

NU13 

S55 

SNC1 

KP1 

NGN 

6 

CRT 

16 



c 

NU23 

CB1 

SNC2 

KP2 

NGT 

7 

CMODE 

17 



c 

G44 

CB2 

SHW1 

KN1 

NMTYPE 

8 

SLN 

18 



c 

G55 

C66 

SHW2 

KN2 

THICK 

9 

SK 

19 



c 

G66 

CB11 

SHQ1 

THETA 

THCK 

10 

RK 

20 



c 

S31 

CB12 

SHQ2 

EPSLON 

1 NSR 

11 

CRT A 

21 




IWPLICIT REAL (A-H,0-Z) 

REAL RCW 

DIMENSION RCM(150000> 

DATA IN5,IN6,IN7,IN8/5,6,5,6/ 

WAXRCW=150000 

IPR=1 

CALL PZERO ( RCW , WAXRCW/ IPR ) 

READ(IN5,2) NSREG,NWT 
WRITE ( IN6 , 9 )NSREG,NMT 

9 FORMAT (/IX, 213, ' NO. OF REGIONS AND NO. OF MATERIALS’) 

2 FORMAT (20 12) 

92 FORMAT (10E1 0.3) 

N1=1+39*NMT*IPR 

N2=N1+IPR 

N3=N2+IPR 

N4=N3+1 

N5=N4+1 

N6=N5+NSREG 

N7=N6+NSREG 

N8=N7+NSREG 

N9=N8 

RCW ( N3 ) =NSREG 
RCM(N4)=NWT 
X=1.D 00 

C READ EPSLON, CURVTR, NSL, NGN, NGT 

CALL METPR0(RCW,IN5,IN6,IN7,IN8,NSREG,NHT,RCM(N5) ,RCW(N6) ,RCW(N7) ,- 
•1,X,1,1,X,IPR,RCW(N1),RCW(N2),1,1,1,1) 

DO 1 I=1,NSREG 

C READ NMTYPE 

CALL HETPR0(RCM,IN5,IN6,IN7,IN8.NSREG,NMT.RCM(N5),RCM(N6>,RCW(N7),- 
•2,X,I,RCW(N9) ,X, IPR, X,X, 1,12, 1,1) 

1 N9=N9+I2 

N10=N9 

DO 3 1=1 , NS REG 
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0006500 

0006600 

0006700 

0006600 

0006900 

0007000 

0007100 

0007200 

0007300 

0007400 

0007500 

0007600 

0007700 

0007800 

0007900 

0008000 

0008100 

0008200 

0008300 

0008400 

0008500 

0008600 

0008700 

0008800 

0008900 

0009000 

0009100 

0009200 

0009300 

0009400 

0009500 

0009600 

0009700 

0009800 

0009900 

0010000 

0010100 

0010200 

0010300 

0010400 

0010500 

0010600 

0010700 

0010800 

0010900 

0011000 

0011100 

0011200 

0011300 

0011400 

0011500 

0011600 

0011700 

0011800 


C READ THICK 

CALL METPR0(RCM,IN5,IN6,IN7,IN8,NSREG,NMT,RCM(N5>,RCn<N6>,RCM<N7>, 
*3,RCtl(N10),I,l,X,IPR,X,X,l,I2,l,l) 

3 N10=N10+I2*IPR 

C READ THCK 

CALL METPRO (RCM, IN5»IN6,IN7,IN8, NSREG, NMT, RCM(N5 ), RCM(N6), RCM(N7 ) , 
*6,RCM(N9),1,1,RCM(N10),IPR,X,X,1,1,1,1> 

N11=N10+NSREG*IPR 
C.....READ MECHANICAL PROPERTIES 

CALL METPR0(RCM,IN5,IN6,IN7,IN8,NSREG,NMT,1,1,1,4,X,1,1,X,IPR,X,X, 
*1,1, 1,1) 

N12=N11+NSREG 

N13=N12+NSREG 

N14=N13+NSREG 

N15=fvl4+NSREG 

N16=N15+NSREG 

N17L1=0 

N17L2=0 

N18L=0 

DO 15 1=1, NSREG 

CALL METPRO (RCM, IN5, IN6, IN7 , IN8 , NSREG , NMT , RCM(N5 ) , RCIKN6 ) , RCM(N7 ) , 
•5,X,1,1,X,IPR,X.X,I»I2,I,I4) 

IF(I4.EQ.0.0R.I4.EQ. T.0R.I4.GT.I ) GOTO 16 
GOTO 15 

16 N17Ll=N17Ll+8*I2-2 
IF(4*I2-2.GT.N57L2) N17L2=4*I2-2 
N18L=N18L+4*I2-2 

15 CONTINUE 

N17=N16+N17L1*2*IPR 

N18=N17+N17L1*N17L2*4*2*IPR 

N19=N18+N18L*IPR 

N20=N19+N17L1*N17L1*IPR 

N21=N20+N17L1*IPR 

N22=N21+N17L1*IPR 

N23=N22+N17L1*IPR 

N24=N23+N17L1*IPR 

N25=N24+N17L1 

URITE(IN6, 17) N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15, 
•N16 , N17 , N17L1 , N17L2 , N18 , N19 , N20 , N21 , N22 , N23 , N24 , N25 

17 FORMAT(20I5) 

CALL HPRG(IN5,IN6,IN7,IN8, NSREG, NUT, RCM(N5),RCM(N6),RCM(N7), 
•RCM(N8),RCM<N9),RCM(N10),RCH(N11),RCM(N12),RCM(N13),RCM(N14), 
•RCn(N15),RCH(N16),RCH(N17),RCtt(N18),RCM(N19>,RCn(N20),RCM(N21), 
•RCM (N22 ) , RCtl ( N23 ) , RCM (N24 ) , N17L1 , N17L2, N18L, N25, MAXRCM, RCH, 

•IPR) 

STOP 

END 

SUBROUTINE METPRO (RCM, IN5,IN6,IN7,IN8, NSREG, N,NSL, NGN, NGT,NXX 
•,THICK,NR,NMTTPE,THCK, IPR, EPSLON,CURVTR, II , 12, 13, 14 ) 

IMPLICIT REAL <A-H,0-2) 

REAL RCM 

DIMENSION RCM(1),NSL(1),NGN(1),NGT(1),THICK(1),NMTYPE(1),THCK(1) 
G0T0(1»2,3,4,5,6),NXX 
1 READ (INS, 17) EPSLON , CURVTR 
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0011900 
0012000 
0012100 
0012200 11 
0012300 
0012400 2 
0012500 
0012600 
0012700 
0012800 
0012900 
0013000 3 
0013100 
0013200 
0013300 
0013400 
0013500 
0013600 6 
0013700 
0013800 
0013900 
0014000 
0014100 
0014200 
0014300 18 
0014400 
0014500 
0014600 
0014700 4 
0014800 
0014900 
0015000 
0015100 
0015200 
0015300 
0015400 
0015500 
0015600 
0015700 
0015800 5 
0015900 
0016000 
0016100 172 
0016200 17 
0016300 14 
0016400 19 
0016500 12 
0016600 13 
0016700 
0016800 16 
0016900 
0017000 
0017100 
0017200 


URITE ( IN6 , 19 ) EPSLON , CURVTR 
DO 11 NR=1 ,NSREG 

READ(IN5,12) NSL(NR) , NGN(NR) ,NGT(NR) 

WRITE (IN6, 13)NR,NSL(NR),NGN(NR),NGT (NR) - 
RETURN 
NL=NSL(NR) 

I2=NL 

WRITE(IN6,14) NR 

READ (INS, 12) (NMTYPE(I),I=1,NL) 

WRITE (IN6, 12) (NMTYPE(I),I=1,NL) 

RETURN 

NL=NSL(NR) 

I2=NL 

WRITE (IN6, 16) NR 
READ(IN5,17) (THICK(I),I=1,NL) 

WRITE (IN6, 17) (THICK (I),I=1,NL) 

RETURN 

J=0 


DO 18 K=1,NSREG 
NL=NSL(K) 

THCK(K)=0.D 00 
DO 18 1=1, NL 
J=J+1 

THCK ( K ) =THCK ( K ) +THI CK ( J ) 

CONTINUE 
WRITE (IN6, 172) 

WRITE (IN6, 17) (THCK(I).I=1,NSREG) 

RETURN 

I=N*IPR 

CALL MPR01 (RCM( 1 ) ,RCM(I+1 ) ,RCM(2*I+1 ) ,RCM(3*I+1 ) ,RCH(4*I+1 ) , 

•RCM (5*1+1 ), RCM (6*1+1 ) ,RCM(7*I+1 ) ,RCM(8*I+1 ) ,RCH(9*I+1 ),RCM( 10*1+1 ) 
*,RCM( 11*1+1 ) ,RCH( 12*1+1 ) ,RCtt( 13*1+1 ) , RCH( 14*1+1 ) ,RCH( 15*1+1) 

• , RCM ( 16*1+1 ) , RCH ( 17*1+1 ) , RCH ( 18*1+1 ) , RCH ( 19*1+1 ) , RCH (20*1+1 ) 
»,RCn(21*I+l ) ,RCH(22*I+1 ) ,RCM(23*I+1 ) ,RCM(24*I+1 ) ,RCH(25*I+1 ) 
*,RCn(26*I+l>,RCn(27*I+l),RCM(28*I+l),RCl1(29*I+l),RCH(30*I+l) 

*, RCM( 31*1+1 ) ,RCtl( 32*1+1 ) ,RCtt( 33*1+1 ) ,RCtt( 34*1+1 ) ,RCM( 35*1+1 ) 
•,RCM(36*I+1 ) ,RCH(37*I+1 ) ,RCM(38*I+1 ) 

*,N,IN5,IN6,IN7,IN8) 

RETURN 
I2=NSL(I1 ) 

I4=NGT(I3) 

RETURN 

FORMAT (/IX, 'SUBLAMINATE THICKNESSES' ) 

FORMAT (8E1 0.4) 

FORMAT (/IX, 'REGION = ',13,', MATERIAL FOR EACH LAYER') 

FORMAT(/1X,2E10.4, ' EPSLON, CURVTR'/) 

FORMAT (2014) 

FORMATdX, 'REGION = ',13,', NO. OF LAYERS = ',13,', NGN = 

•' ,13, ' , NGT = ',13) 

FORMAK/1X, 'REGION = ',13,', PLY THICKNESSES') 

END 

SUBROUTINE MPRG ( IN5 , IN6 , IN7 , INS , NSREG, NMT , NSL , NGN , NGT , 

•NMTYPE , THICK, THCK, NSR, NZR, NRR, NCR , 

•NBCR , CRT, CMODE , SLN , SK, RK, CRTA, 
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0017300 

00174100 

0017500 

0017600 

0017700 

0017800 

0017900 

0018000 

0018100 

0018200 

0018300 

0018000 

0018500 

0018600 

0018700 

0018800 

0018900 

0019000 

0019100 

0019200 

0019300 

00194100 

0019500 

0019600 

0019700 

0019800 

0019900 

0020000 

0020100 

0020200 

0020300 

00204100 

0020500 

0020600 

0020700 

0020800 

0020900 

0021000 

0021100 

0021200 

0021300 

0021400 

0021500 

0021600 

0021700 

0021800 

0021900 

0022000 

0022100 

0022200 

0022300 

0022400 

0022500 

0022600 


•CRTS, CRTC , I SGN , NCN1 , NCN2 , NOW , N25 , MAXRCM . RCH ,IPR) 

IMPLICIT REAL (A-H,0-Z> 

REAL RCM 

DIMENSION NSL(1>,NGN(1),NGT(1),NZR(1),NRR(1>,NCRU>,NBCR(1>,NSR<1> 
•,CRT(NCN1,1) ,CM0DE(NCNl,NCN2,l) r SLN(l),CRTA(l),THCK(l), 

•CRTB ( 1 ) , CRTC ( 1 ) , SK ( NCN1 , 1 ) , RK ( 1 ) , RCM ( 1 ) , I SGN ( 1 ) » THICK ( 1 ) , NMTYPE ( 1 ) 
CALL PZER0(CRT,NCN1*2) 

CALL PZERO ( CMODE , NCN1 *NCN2*4 ) 

CALL PZERO (SLN,NCN4> 

NG1=0 

NG2=0 

NG4=0 

DO 1244 NR=1,NSREG 
NRGN=NGN(NR) 

NZR(NR)=0 

NRR(NR)=0 

NCR(NR)=0 

NSR(NR)=0 

NL=NSL(NR) 

H=THCK(NR> 

NPLY=NL 

IF(NR.GT.l) NG1=NG1+NSL(NR-1) 

NRT=0 

NCT=0 

IF (NGT(NR) .EQ.0) GOTO 12345 
IF(NGT(NR).EQ.NR) GOTO 12345 
NR1=NGT(NR) 

NZR(NR)=NZR(NR1) 

NSR(NR)=NSR(NR1) 

NRR(NR)=NRR(NR1) 

NCR(NR)=NCR(NR1 ) 

NBCR ( NR ) =NBCR ( NR1 ) 

GOTO 1244 

12345 DO 1245 NSYM=1,3 
NSTREI=5*NPLY-1 
N25L=4*NSL(NR)-2 
N26=N25+N25L**2*IPR 
N27 =N26+N25L»*2* IPR 
N28=N27+N25L**2* IPR 
N29=N28+N25L**2*IPR 
N30=N29+N25L**2*IPR 
N31=N30+NSTREI*NSTREI*IPR 
N32=N31+NSTREI*IPR 
N33=N32+ (4*NL-2)*IPR 
N34=N33+(4*NL-2>*1 

IF < N34 . GT . MAXRCM ) WRITE (IN6, 39846 ) N34 , MAXRCM 
IF(N34.GT.MAXRCH> STOP 

CALL DEQHAT(NSYM,NL,NEQ,NBC r NGl,NG2,NR,NMT,NSTREI,H, 

•RCM ( 1+21 »NMT ) r RCM ( 1+31«NMT ) , RCM ( 1+33*NMT ) , RCM ( 1+35*NMT ) , 

•RCM C 1 +37*NMT ) , THICK , SLN, CRT , CMODE , RCM ( N25 ) , 

•RCM(N26) ,RCM(N27),RCM(N28) ,RCH(N29),RCM(N30),RCM(N31),RCM(N32) , 
•RCM ( N33 ) , RCM , 

•NCN1 , NCN2 , NSREG , NMTYPE , NG4 , N25L , IPR ) 

IF(NSYM.GT.l) GOTO 2356 
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0022700 

0022800 

0022900 

0023000 

0023100 

0023200 

0023300 

0023400 

0023500 

0023600 

0023700 

0023800 

0023900 

0024000 

0024100 

0024200 

0024300 

0024400 

0024500 

0024600 

0024700 

0024800 

0024900 

0025000 

0025100 

0025200 

0025300 

0025400 

0025500 

0025600 

0025700 

0025800 

0025900 

0026000 

0026100 

0026200 

0026300 

0026400 

0026500 

0026600 

0026700 

0026800 

0026900 

0027000 

0027100 

0027200 

0027300 

0027400 

0027500 

0027600 

0027700 

0027800 

0027900 

0023000 


NXEUQJ=4*NSL(NR)-2 
WRITE (IN6, 1459) NR 

WRITE ( IN6 , 1458 ) ( SLN (NG4+I ) , 1=1 , NXEUQJ ) 

IF (NRGN.EQ.O )NBCR(NR)=(8*NSL(NR)-2) 

IF ( NRGN . NE . 0 ) NBCR ( NR ) = ( 8*NSL ( NR ) -2-6 ) /2+6 

NSR(NR)=NXEUQJ 

NG4=NG4+1 

I1=NG2+1 

NZT=6 

I2=NG2+NZT 

DO 191 1=11,12 

CRT ( I , 1 ) =0 . D 00 

CRT(I,2)=0.D 00 

WRITE (IN6, 1463) I, NR 

DO 191 J=l, NXEUQJ 

191 WRITE (IN6, 870) (CMODE(I, J,K),K=1,4) 

NZR ( NR ) =NZR ( NR )+NZT 
NG2=NG2+NZT 

2356 IF(NPLY.EQ.l) GOTO 1245 
DO 1369 NIRT=1,2 
WRITE ( IN8 , 2357 )NR, NSYN 
READ ( IN7 , 2 ) NK10 
IF(NIRT.EQ.l) WRITE (IN8, 1367) 

IF(NIRT.EQ.2) WRITE (IN8, 1370) 

READ ( IN7 , 2 ) J 
IF(J.NE.O) GTT9 1369 
IF(NIRT.NE.l) GOTO 4884 
N25L=4*NSL ( NR ) -2 
N26 =N25+N25L**2* I PR 
N27 =N26+N25L**2* I PR 
N28=N27+N25L**2*IPR 
N29=N28+N25L**2* IPR 
N30=N29+N2SL**2*IPR 
N31=N30+N25L *IPR 
N32=N31+N25L*IPR 
N33=N32+N25L**2*IPR 

IF (N33 . GT . nAXRCM ) WRITE ( IN6 , 39845 )N33 , HAXRCH 
IF(N33.GT.HAXRCM) STOP 

39845 FORMAT ( IX, 'N33 = ',17,', WAXRCM = »,I7> 

CALL RROOT(NSYM,NEQ,NT,NG2,NRGN,NK10,NRT,NCT,IN7,IN8, 
•CRT,CH0DE,RCM(N25) ,RCtl(N26) ,RCM(N27) ,RCM(N28),RCn(N29) # CRTA, 
•RCM ( N30 ) , RCM < N31 ) , RCM ( N32 ) , NCN1 , NCN2 , N25L ) 

GOTO 4885 

4884 CONTINUE 

N25L=4*NSL ( NR ) -2 

N26=N25+N25L**2*IPR 

N27=N26+N25L**2*IPR 

N28=N27+N25L*»2*IPR 

N29=N28+N25L**2*IPR 

N30 =N29+N25L**2*I PR 

N31=N30+N25L*2*IPR 

N32=N31+N25L«2»IPR 

N33=N32+N25L**2*2*IPR 

N34=N33+N25L*IPR 
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0028100 

0028200 

0028300 

0028400 

0028500 

0028600 

0028700 

0028800 

0028900 

0029000 

0029100 

0029200 

0029300 

0029400 

0029500 

0029600 

0029700 

0029800 

0029900 

0030000 

0030100 

0030200 

0030300 

0030400 

0030500 

0030600 

0030700 

0030800 

0030900 

0031000 

0031100 

0031200 

0031300 

0031400 

0031500 

0031600 

0031700 

0031800 

0031900 

0032000 

0032100 

0032200 

0032300 

0032400 

0032500 

0032600 

0032700 

0032800 

0032900 

0033000 

0033100 

0033200 

0033300 

0033400 


IF (N34 .GT. MAXRCM ) WRITE ( IN6 , 39846 )N34 , MAXRCM 
39846 FORMAT ( IX, 'N34 = ,»,I7,», MAXRCM = »,I7) 

IF (N34.GT. MAXRCM) STOP 

CALL CROOT(NSYM,NEQ,NT,NG2,NRGN,NK10,NRT,NCT,IN7,IN8, 
•CRT , CMODE , RCM ( N25 ) , RCtt ( N26 ) , RCM ( N27 ) , RCM ( N28 ) , RCM ( N29 > , CRTA, 
•RCM(N30 ) ,RCH(N31 ) ,RCM(N32 ),RCM(N33),NCN1,NCN2,N25L) 

4885 IF (NSYH.EQ.3) GOTO 3241 

IF(NIRT.EQ.l.AND.NRT.LT.l) GOTO 3241 
IF (NIRT.EQ.2. AND. id .LT.2) GOTO 3241 
NEVN=(HPLY+1 >/2 
N2EVN=(4*NPLY-2)/2+2 
N3EVN= (4*NPLY-2 )/2+l 
N4EVN=(4*NPLY-2) 

IF (NPLY/2*2 . NE . NPLY . AND . NSYM. EQ . 1 ) I3=N3EVN-1 
IF (NPLY/2*2 . NE . NPLY . AND . NSYM. EQ . 2 ) I3=N3EVN 
IF(NPLY/2«2.EQ.NPLY .AND. NSYM. EQ.l) I3=N3EVN 
IF (NPLY/2*2 . EQ . NPLY . AND . NSYM . EQ . 2 ) I3=N3EVN- 1 
DO 7119 I=1,NEVN 
J=4* (NPLY+l-I ) 

ISGN(J-3)= <4*I-3)*(-l.)»*CNSYM) 

ISGN(J-2>= (4*1-2 )*C-1.)**(NSYM+1) 

ISGN(«J-5)= (4*1-1 )* ( -1 . )•* (NSYM+1 ) 

7119 ISGN(J-4)= (4*1-0 )*(-l.)**(NSYM> 

I1=NG2+1 

IF (NIRT.EQ.l) I2=NG2+NRT 
IF(NIRT.EQ.2) I2=NG2+NCT 
DO 424 1=11,12 
DO 424 L =1,4 
DO 7126 J=N2EVN , N4EVN 
K=ISGN(J) 

JSGN=K/IABS(K) 

K=IABS(K) 

7126 CMODE(I,J,L)=JSGN*CMODE(I,K,L> 

424 CMODE(I,I3,L)=O.D 00 

3241 IF(NIRT.EQ.l) WRITE (IN8, 1323) 

IF (NIRT.EQ.2) WRITE (IN8, 1324) 

READ (IN7,2 )J 
IF(J.NE.O) GOTO 1369 
IF(NIRT.EQ.l.AND.NRT.LT.l) GOTO 1369 
IF(NIRT.EQ. 2. AND.NCT.LT. 2) GOTO 1369 
I1=NG2+1 

IF (NIRT.EQ.l) I2=NG2+NRT 
IF(NIRT.EQ.2) I2=NG2+NCT 
NXEUQJ=4*NSL(NR)-2 
DO 880 13=11,12 

WRITE (IN6, 8610) I3,CRT(I3,1),CRT(I3,2) ,NR 
DO 880 <J=1,NXEWQJ 

880 WRITE (IN6, 870) (CMODE(I3, J,K),K=1,4) 

IF(NIRT.EQ.l) NG2=NG2+NRT 

IF (NIRT.EQ.2) NG2=NG2+NCT 

IF (NIRT.EQ.l) NRR(NR)=NRR(NR)+NRT 

IF (NIRT.EQ.2) NCR(NR)=NCR(NR)*NCT 

NRT=0 

NCT=0 
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0033500 1369 

CONTINUE 

0033600 

WRITE C IN8 , 483 ) 

0033700 

READ(IN7,2)J 

0033800 

IF(J.EQ.O) GOTO 1245 

0033900 

NRTNCT =NRT +NCT 

0034000 

11=1 

0034100 

12=1 

0034200 

13=1 

0034300 6736 

IF (CRT(NG2+I1,2) .NE.O. ) GOTO 6734 

0034400 

CRTA ( 12+0 ) =CRT ( NG2+I 1,1) 

0034500 

CRTA ( 12+1 ) =CRT CNG2+I1 , 1 ) 

0034600 

11=11+2 

0034700 

12=12+2 

0034800 

GOTO 6735 

0034900 6734 

DO 6847 1=1,4 

0035000 

CRTB( I3+I-1 ) =CRT (NG2+I1 , 1 ) 

0035100 6847 

CRTC ( I 3+ I - 1 ) =CRT ( NG2+ I 1 , 2 ) 

0035200 

13=13+4 

0035300 

11=11+4 

0035400 6735 

IF ( 11. GT. NRTNCT) GOTO 6737 

0035500 

GOTO 6736 

0035600 6737 

IF(NRT.LT.2) GOTO 6589 

0035700 

DO 6738 11=1, NRT 

0035800 

CRT ( NG2+I 1 , 1 ) =CRTA ( I 1 ) 

0035900 6738 

CRT CNG2+I1 , 2 )=0 . D 00 

0036000 6589 

IF (NCT.LT.4) GOTO 6588 

0036100 

DO 6739 11=1, NCT 

0036200 

CRT (NG2+NRT+I1 , 1 ) =CRTB( 11 ) 

0036300 6739 

CRT ( NG2+NRT+1 1,2) =CRTC (ID 

0036400 6588 

CONTINUE 

0036500 

DO 1906 1=1, NRT 

0036600 1906 

WRITE (IN8, 1908) NRT,I,CRT(NG2+I,1),CRT(NG2+I,2) 

0036700 

DO 1907 1=1, NCT 

0036800 1907 

WRITE ( IN8 , 1908 ) NCT , I , CRT ( NG2+NRT+1 , 1 ) , CRT (NG2+NRT+1 , 2 ) 

0036900 

WRITE (IN8, 1943) 

0037000 

READ ( IN7 , 2 )NRTXL 

0037100 

IF(NRTXL.LT.O) NRT=0 

0037200 

IF(NRTXL.GT.O) NRT =NRTXL 

0037300 

WRITE (IN8, 1944) 

0037400 

READ ( IN7 , 2 ) NCTXL 

0037500 

IF(NCTXL.LT.O) NCT=0 

0037600 

IF (NCTXL. GT.O) NCT=NCTXL 

0037700 

WRITE(IN8,1909) 

0037800 

READ(IN7,2) 0 

0037900 

IF(J.EQ.O) GOTO 2356 

0038000 

IF (NSYtl. EQ . 1 ) NT=( (NBC-1 ) -NCT-NRT-4+2 ) /2 

0038100 

IF(NSYM.EQ.2) NT=( <NBC-l)-NCT-NRT-2+2)/2 

0038200 

IF (NSYtl. EQ. 3) NT=( (NBC-0 )-NCT-NRT-6+2)/2 

0038300 

CALL PZER0(CRTA,NCN1) 

0038400 

NK10=0 

0038500 

N25L=4»NSL(NR)-2 

0038600 

N26=N25+N25L**2*IPR 

0038700 

N27 =N26+N25L**2*IPR 

0038800 

N28=N27+N25L**2*IPR 
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0038900 N29=N28+N25L**2*IPR 

0039000 N30=N29+N25L**2*IPR 

0039100 N31=N30+N25L*<*2*IPR 

0039200 CALL COEFF1 ( NSYH, CRTA, NEQ , NCT , NRT , NRGN , NT , NG2 , RCM ( N31 ) , RCtl ( N25 ) , - 

0039300 •RCM(N26),RCM(N27),RCH(N28),RCn(N29),CRT r NK10,NCNl,RCM<N30),N25L) 

0039400 HRITE (IN8, 44844 )( (I,CRTA(I)),I=1,NT) 

0039500 GOTO 2356 

0039600 1245 CONTINUE 
0039700 1244 CONTINUE 

0039800 HRITE (IN6, 1654) (NSL(I>,I=1,NSREG) 

0039900 HRITE (IN6, 1656) (NZR(I>,I=1,NSREG) 

0040000 HRITE (IN6, 1657) (NRR(I),I=1,NSREG) 

0040100 HRITE (IN6, 1658) (NCR(I),I=1,NSREG) 

0040200 URITE(IN6,1659) (NBCR(I),I=1,NSREG> 

0040300 HRITE (IN6, 1652) <NSR(I),I=1,NSREG) 

0040400 NTBC=0 

0040500 DO 1411 I=1,NSREG 

0040600 1411 NTBC=NTBC+NZR(I )+NRR(I >+NCR(I ) 

0040700' HRITE (IN6, 8603) NTBC 

0040800 CALL PZERO (SK, NCN1*NCN1 ) 

0040900 CALL PZERO i RK, NCN1 ) 

0041000 DO 1413 m=l,NTBC 

0041100 IF(ni.EQ.l) HRITE (IN6, 8606) 

0041200 112=1 

0041300 14442 READ(IN5,8605)NBCREG,NBCPLY,NBCV,YNBC,ZNBC,PNBC,NBCRG1 
0041400 C IF (H2.EQ.) > HRITE ( I N6, 6667) 

0041500 HRITE C IN6 , 8605 ) NBCREG , NBCPLY , NBCV , YNBC, ZNBC , PNBC , NBCRG1 

0041600 NR1=NBCREG 

0041700 NR=NGT (NR1 ) 

0041800 IF(NR.EQ.O) NR=NR1 

0041900 NRGN=NGN(NR) 

0042000 NPLY=NBCPLY 

0042100 NV=NBCV 

0042200 PS = PNBC 

0042300 SSSGN=1. 

0042400 SSSGN1=1. 

0042500 IF012.GT.1) SSSGN=PS 

0042600 IF(M2.GT.l) SSSGN1=PS 

0042700 IF(SSSGN.EQ.O.) SSSGN=1. 

0042800 IF(SSSGN1.EQ.0.) SSSGN1=1. 

0042900 IF(W2.GT.l) PS=0.D 00 

0043000 NV=IABS(NV) 

0043100 Y=YNBC 

0043200 Z=ZNBC 

0043300 H=THCK(NR) 

0043400 NL=NSL(NR) 

0043500 NV1=NV 

0043600 IF(NRGN.NE.l) GOTO 64851 

0043700 IF(DABS(Y).LT.1.E-10) GOTO 64851 

0043800 IF(NVl.NE.l) GOTO 64852 

0043900 NPLY=1 

0044000 GOTO 64851 

0044100 64852 IF(NV1.NE.3) GOTO 64854 

0044200 NPLY=1 
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0044300 

0044400 

0044500 

0044600 

0044700 

0044800 

0044900 

0045000 

0045100 

0045200 

0045300 

0045400 

0045500 

0045600 

0045700 

0045800 

0045900 

0046000 

0046100' 

0046200 

0046300 

0046400 

0046500 

0046600 

0046700 

0046800 

0046900 

0047000 

0047100 

0047200 

0047300 

0047400 

0047500 

0047600 

0047700 

0047800 

0047900 

0048000 

0048100 

0048200 

0048300 

0048400 

0048500 

0048600 

0048700 

0048800 

0048900 

0049000 

0049100 

0049200 

0049300 

0049400 

0049500 

0049600 


NV=1 

GOTO 64851 

64853 NPLY=1 
NV =3 

GOTO 64851 

64854 IF (NV1 .NE.5) GOTO 64851 
NV=5 

NPLY=1 
64851 NG1=0 
NG2=0 
NG3=0 
NG4=0 

NCT=NCR(NR) 

NZT=NZR(NR) 

KRT=NRR(NR) 

NST=NSR(NR> 

DO 1414 1=1 ,NR1 
IFCI.EQ.l) GOTO 1414 
IF(I.GT.NR) GOTO 1414 
NG1=NG1+NSL(I-1> 

IF (NGT(I ) .EQ.O .OR.NGTd ) .EQ.NR)NG2=NG2+NZR(I - 1 )+NRR(I-l)+NCRd _ l ) 
IF ( NGT ( I > . EQ . 0 . OR . NGT ( I ) . EQ . NR )NG4=NG4+NSR ( I - 1 ) 

NG3=NG3+NBCR(I-1> 

1414 CONTINUE 

NG1=NG1+NPLY 

IF(NRGN.NE.l) GOTO 66710 

IF(DABS(Y>.LT.1.E-10> GOTO 66710 

IF (NV1 .NE.3) GOTO 66710 

IF (NV.EQ.3) SSSGN =SSSGN1 

IF (NV .NE.l) GOTO 66710 

ECCN=H/2 

DO 66711 ICFG54=1,NPLY 
66711 ECCN=ECCN-THICK(NG1-NPLY+ICFG54) 

ECCN=ECCN+THICK (NG1 )/2 
SSSGN=SSSGN1*ECCN 
66710 I111=NG3+1 

I112=NG3+NBCR(NR) 

DO 1415 N1=I111,I112 
NG2=NG2+1 

IF(CRT(NG2,1).EQ.0..AND.CRT(NG2,2).EQ.0.) HT=1 
IF(CRT(NG2,l).NE.O..ftND.CRT(NG2,2).EQ.O.) HT=2 
IF (CRT(NG2, 1) .NE.O. .RND.CRT(NG2,2) -NE.O- ) HT=3 
NTM=IPR*NMT 

CALL FT)N(NL,NPLY,NG1,NG2,NG4,NRGN,NV,HT,N11T,H,Y,Z,F,FPS, 

•RCll ( 1 + 1 8*NTN ) , RCM (1+1 9*NTM ) , RCM ( 1 +21 *NTM ) , RCW ( 1 +25*NT11 ) , 

•RCM( 1 +27*NT!1 ) , RCW ( 1+29*NT«) , RCM ( 1+31*NTM ) , RCM ( 1+33*NTM) , 

•RCM ( 1+35*NTM ) , RCM ( 1+37*NTM ) , RCM ( 1+39*NTM ) , RCM <2+39*NTM ) , 

•THICK ,NMTYPE, 

•SLN , CRT , CMODE , RCM, NCN1 , NCN2 , IPR ) 

FPS1=FPS 

SK ( m , N1 ) =SK (Ml , N1 ) +SSSGN*F 

C WRITE (IN6, 6666) NR,NPLY,NG1,NG2,MT,M1,N1,NV,Y,Z,F,FPS»SSSGN,PS 

6667 FORMAT (/IX, 'NR,NPLY,NG1,NG2,MT.M1,N1,NV,Y,Z,F,FPS, SSSGN, PS' ) 

6666 FORMAT (IX, 813, 6E9.-J) 
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0049700 

0049800 

0049900 

00S0000 

0050100 

0050200 

0050300 

0050400 

0050500 

0050600 

0050700 

0050800 

0050900 

0051000 

0051100 

0051200 

0051300 

0051400 

0051500 

0051600 

0051700 

0051800 

0051900 

0052000 

0052100 

0052200 

0052300 

0052400 

0052500 

0052600 

0052700 

0052800 

0052900 

0053000 

0053100 

0053200 

0053300 

0053400 

0053500 

0053600 

0053700 

0053800 

0053900 

0054000 

0054100 

0054200 

0054300 

0054400 

0054500 

0054600 

0054700 

0054800 

0054900 

0055000 


1415 CONTINUE 

IF(NRGN.NE.l) GOTO 14140 
IF(DABS(Y).LT.1.E-10> GOTO 14140 
IFCNV1.NE.1 ) GOTO 69941 
NPLY=NPLY+1 

RK(M1 >=RK(M1 >-SSSGNoFPSl 

IF(M2.EQ.l. AND.NPLY.EQ.NL) RK(M1 >=RK(M1 )-PS 
IF (NPLY.GT.NL) GOTO 14131 
GOTO 64851 

69941 IF (NV1 .NE.3) GOTO 69942 
NPLY=NPLY+1 
IF(NV.NE.l) GOTO 69943 
RK(M1)=RK(H1)-SSSGN*FPS1 

IF(M2.EQ.l. AND.NPLY.EQ.NL) RK(M1 >=RK(M1 )-PS 
IF (NPLY.GT.NL) GOTO 64853 
GOTO 64851 

69943 RK(M1)=RK(M1)-$SSGN*FPS1 

IF(H2.EQ.l. AND.NPLY.EQ.NL) RK(M1 >=RK(M1 )-PS 
IF(NPLY.GT.NL) GOTO 14131 
GOTO 64851 

69942 IF (NV1 .NE.5) GOTO 14140 
NPLY=NPLY+1 

RK (Ml > =RK (Ml ) -SSSGN*FPS1 

IF(M2.EQ.l. AND.NPLY.EQ.NL) RK(M1 )=RK(M1 )-PS 
IF (NPLY.GT.NL) GOTO 14131 
GOTO 64851 

14140 RK(M1)=RK(M1)-SSSGN*FPS1 

IF (M2.EQ.1 ) RK(M1)=RK(M1)-PS 
14131 M2=M2+1 

IF(NBCRGl.EQ.O) GOTO 1413 
GOTO 14442 
1413 CONTINUE 

CALL SKINV ( SK, NCN1 , SK , NCN1 , RK , NCN1 , 1 , SIMUL , NTBC > 

386 READ ( I N5 , 4 02 ) NDREG , NDPLY , NDV , NDVT , YND , ZND , STPND , NPTS 
IF(NDREG.EQ.O) GOTO 387 
WRITE (IN6, 339) 

Y=YND 

Z=2ND 

N=1 

704 FD=0.D 00 
NR1=NDREG 
NR=NGT(NR1) 

IF (NR.EQ.O) NR=NR1 
NRGN=NGN(NR) 

NPLY=NDPLY 

NV=NDV 

H=THCK(NR) 

NL=NSL(NR) 

NG1=0 

NG2=0 

NG3=0 

NG4=0 

NCT=NCR(NR) 

NZT=NZR(NR) 
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0055100 NRT=NRR(NR) 

0055200 NST=NSR(NR) 

0055300 DO 373 J=1,NR1 

0055400 IF(J.EQ.l) GOTO 373 

0055500 IF (J.GT.NR) GOTO 373 

0055600 NG1=NG1+NSL(0-1) 

0055700 IF(NGT(J).EQ.O.OR.NGT(J).EQ.NR)NG2=NG2+NZR(J-1)+NRR(J-1)+NCR(J-1) 

0055800 IF(NGT(J).EQ.O.OR.NGT(J).EQ.NR)NG4=NG4+NSR(J-1) 

0055900 NG3=NG3+NBCR(J-1) 

0056000 373 CONTINUE 
0056100 NG1=NG1+NPLY 

0056200 I111=NG3+1 

0056300 I112=NG3+NBCR(NR> 

0056-300 DO 703 N1=I111,I112 

0056500 NG2=NG2+1 

0056600 IF(CRT(NG2,1).EQ.0..AND.CRT(NG2,2).EQ.0.) MT=1 

0056700 IF(CRT(NG2,1).NE.0..AND.CRT(NG2,2).EQ.0.) MT=2 

0056800 IF(CRT(NG2,l).NE.O..AND.CRT(NG2,2).NE.O.) MT=3 

0056900 NTM=IPR*NMT 

0057000 CALL FUN(NL,NPLY,NG1,NG2,NG4,NRGN,NV,MT,NMT,H, Y,Z,F,FPS, 

0057100 »RCM ( 1+16+NTM ) , RCH ( 1+19*NTM > , RCM ( 1+21*NTM ) , RCM ( 1+25*NTM ) , 

0057200 *RCM(1+?7*NTM),RCM(1+29*NTM),RCM(1+31*NTM),RCM(1+33*NTM), 

0057300 •RCM(1+35*NTM),RCM(1+37*NTH),RCM(1+39*NTM),RCM(2+39*NTM), 

0057400 »THICK ,N«TYPE, 

0057500 «SLN , CRT , CMODE , RCM , NCN 1 , NCN2 , I PR ) 

0057600 FPS1=FF~ 

0057700 IF(N1.EQ.I112) FD=FD+FPS1 

0057800 703 FD=FD+RK(N1 )*F 

0057900 HRITE(IN6,338) Y,FD,Z,NR,NPLY,NV 

0058000 N=N+1 

0058100 IF (N.GT. 1.AND.NDVT.EQ.2) Y=Y+STPND 

0058200 IF(N.GT.1.A!®.NDVT.EQ.3) Z=Z+STPND 

0058300 IF(N.LE.NPTS+1) GOTO 704 

0058400 GOTO 386 

0058500 387 CONTINUE 

0058600 1652 F0RMAT('1',/1X, ' NO. OF VARIABLES IN EACH REGION r ,8(lX,I3>) 

0058700 1654 FORMAT (IX , ' NO. OF PLYS IN EACH REGION r ,8(lX,I3>) 

0058800 1656 FORMAK1X, * NO. OF ZERO ROOTS IN EACH REGION r ,8(lX,I3>) 

0058900 1657 FORMAT (IX , ' NO. OF REAL ROOTS IN EACH REGION »,8(1X,I3)) 

0059000 1658 FORMAT (IX, r NO. OF COMP. ROOTS IN EACH REGION r ,8(lX,I3)) 

0059100 1659 FORMAK1X, r NO. OF B. CONDITIONS IN EACH REGION' ,8(1X, 13) ) 

0059200 1909 FORMAT (IX, 'GIVE 0 IF POLYNOMIAL IS NOT DESIRED') 

0059300 1943 FORMAT (IX, 'GIVE NRT FOR POLYNOMIAL. .. .DEFAULT IS AS IS'/IX, - 

0059400 *'A NEGATIVE VALUE INITIATES NRT TO O') 

0059500 1944 FORMAT (IX, 'GIVE NCT FOR POLYNOMIAL DEFAULT IS AS IS'/IX, - 

0059600 •' A NEGATIVE VALUE INITIATES NCT TO O') 

0059700 339 F0RMAT(/1X, 'Y,F,Z,NR,NPLY,NV' ) 

0059800 8606 F0RMAT('1',/1X, 'NBCREG,NBCPLY,NBCV,YNBC,ZNBC,PNBC,NBCREG' ) 

0059900 8603 FORMAT (/IX, 14, ' = NTBC' ) 

0060000 483 F0RMAT(/1X, 'GIVE 0 IF NO POLYNOMIAL DESIRED') 

0060100 1324 FORMAT (IX, 'GIVE 0 IF COMPLEX ROOTS ARE TO BE STORED') 

0060200 1370 FORMAT (IX, 'GIVE 0 IF YOU UANT COMPLEX ROOTS') 

0060300 8610 F0RMAT(/1X,I4,1X,2D15.8, ' EIGEN VALUE AND MODE FOR REGION = ',13) 

0060400 870 F0RMAT(6(1X,F17.6) ) 
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0060S00 

0060600 

0060700 

0060800 

0060900 

0061000 

0061100 

0061200 

0061300 

0061400 

0061500 

0061600 

0061700 

0061800 

0061900 

0062000 

0062100 

0062200 

0062300, 

0062400 

0062500 

0062600 

0062700 

0062800 

0062900 

0063000 

0063100 

0063200 

0063300 

0063400 

0063500 

0063600 

0063700 

0063800 

0063900 

0064000 

0064100 

0064200 

0064300 

0064400 

0064500 

0064600 

0064700 

0064800 

0064900 

0065000 

0065100 

0065200 

6065300 

0065400 

0065500 

0065600 

0065700 

0065800 


1323 FORMAT (IX, ’GIVE 0 IF REAL ROOTS ARE TO BE STORED') 

1367 FORMAT (IX, 'GIVE 0 IF YOU UANT REAL ROOTS') 

2357 FORMAT (IX, 'GIVE 0 IF YOU DO NOT UANT POLYNOMIAL APPROX.'/ 

•IX, 'NR = ',13, ' NSYM = ',13) - 

1463 FORMAT (/IX, 14,' ZERO MODE, ’,13,' REGION') 

1459 F0RMAT('1’,/1X,' PARTICULAR SOLUTION FOR REGION = ’,12) 

8605 FORMAT (3I4,3D12.6,3I4,2D12.6) 

402 FORMAT (414, 3D12. 6, 14, 14) 

338 F0RMAT(3(1X,F17.6),3(1X,I3,1X>) 

1908 F0RMAT(1X,2(I3, IX) ,2(D15.8, IX) ) 

44844 FORMAT ( 13, 2X,D23. 16) 

1458 FORMAT ( 10 (1X,D10.4) ) 

123 FORMAT (8(1X,D9.3)) 

12 FORMAT(6(1X,E10.4) ,/3(lX,E10.4) ) 

435 FORMAT (8D1 0.3) 

2 FORMAT (2012) 

RETURN 

END 

SUBROUTINE MPR01(E11,E22,E33,NU12,NU13,NU23,G44,G55,G66,S31,S32, 
•S33 , S44 , S45 , S55 , CB1 , CB2 , CB6 , CB1 1 , CB1 2 , CB16 , CB22 , CB26 , CB66 , SNC1 , 
•SNC2 , SHM1 , SHM2 , SHQ1 , SHQ2 , KQ1 , KQ2 , KM1 , KM2 , KP1 , KP2 , KN1 , KN2 . THETA, 

•N, IN5,IN6, IN7, IN8) 

IMPLICIT REAL (A-H,0-Z> 

REAL KM1,KM2,KN1,KN2,KQ1,KQ2,KP1,KP2,NU12,NU13,NU23 
DIMENSION E11(1),E22(1),E33(1),NU12(1),NU13(1),NU23(1),G44(1), 
•G55(1),G66(1),S31(1),S32(1),S33(1),S44(1),S45(1),S55(1),CP1(1), 
•CB2(1),CB6(1),CB11(1),CB12(1),CB16(1),CB22(1),CB26(1),CB66U), 
•SNC1 (1 ),SNC2(1 ),SHM1 (1 ) ,SHM2(1 ) ,SHQ1 (1 ) ,SHQ2(1 ) ,KQ1 (1 ) ,KQ2(1 ) , 

•KM1 ( 1 ) , KM2 ( 1 ) , KP1 ( 1 ) , KP2 ( 1 ) , KN 1 ( 1 ) , KN2 ( 1 > , THETA ( 1 ) 
•,T1(6,6),T2(6,6),T1I(6,6),ST(6,6),SF (6,6),VC6(6) 

C LAYER FLEXIBILITI MATRICES 

PI=3. 141592653589793 
URITE(IN6, 149) 

149 FORMAT ( ' 1 ' ) 

URITE(IN6,5) 

DO 1242 K=1,N 

READ(IN5,4) E11(K),E22(K),E33(K),NU12(K),NU13(K),NU23(K),G44(K>, 
•G55(K),G66(K),THETA(K) 

1242 URITE(IN6,123)E11(K),E22(K),E33(K),NU12(K),NU13(K),NU23<K),G44(K), 

•G55 ( K ), G66 ( K ), THETA ( K ) 

5 F0RMAT(/1X, 'Ell, E22,E33,NU12,NU13,NU23,G44,G55,G66, THETA’) 

4 FORMAT (10F8. 4) 

DO 101 K=1,N 
CALL PZER0(ST,36) 

ST(1,1)=1./E11(K) 

ST(1,2)=1./E11(K)»NU12(K)»(-1) 

ST(1,3)=1./E11 (K)»NU13(K)»(-1) 

ST(2,1 )=ST(1,2) 

ST(2,2)=1./E22(K) 

ST(2,3)=1./E22(K)*NU23(K)*(-1) 

ST(3,1)=ST(1,3) 

ST(3,2)=ST(2,3) 

ST(3,3)=1./E33(K) 

ST(4,4)=1./G44(K) 
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0065900 


ST(5,5)=1./G55(K) 

0066000 


ST(6,6)=1./G66(K) 

0066100 


CALL PZERO (Tl, 36) 

0066200 


CO=DC0S(THETA(K)*PI/180) 

0066300 


SI=DSIN(THETA(K)*PI/180) 

0066400 


Tl ( 1 , 1 ) =CO*CO 

0066500 


T1(1,2)=SI*SI 

0066600 


Tl(l,6)=-2.*CO*SI 

0066700 


T1(2,1)=SI*SI 

0066800 


T1(2,2)=C0*C0 

0066900 


T1(2,6)=2*C0*SI 

0067000 


Tl(3,3)=l 

0067100 


T1(4,4)=C0 

0067200 


T1(4,5)=SI 

0067300 


Tl (5,4 >=-SI 

0067400 


Tl (5,5)=C0 

0067500 


Tl (6,1 )=CO*SI 

0067600 


Tl (6,2)=-CO*SI 

0067700 


Tl ( 6 , 6 ) =CO*CO-SI*SI 

0067800 


CALL PZERO (T2, 36) 

0067900 


CALL PZERO (SF, 36) 

0068000 


CALL SKINV(T1,6,T2,6,VC6,6,-1,$IMUL,6) 

0068100 


CALL BID0T(ST,6,T2,6,SF,6) 

0068200 


DO 102 1=1,6 

0068300 


DO 102 J=1 ,6 

0068400 

102 

T2(I,J)=T1(I,J) 

0068500 


T2(6,l )=2.*CO*SI 

0068600 


T2(6,2 )=-2.*C0*SI 

0068700 


T2(l,6)=-CO*SI 

0068800 


T2(2,6)=C0*SI 

0068900 


CALL BID0T(T2,6,SF,6,ST,6) 

0069000 


WRITE (IN6, 105)K 

0069100 

105 

FORMAK/1X, 'FLEXIBILITY MATRIX S FOR ’,12/ MATERIAL') 

0069200 


DO 1010 1=1,6 

0069300 

1010 

HRITE(IN6,123)ST(I,1),ST(I,2),ST(I,3),ST(I,4),ST(I,5),ST(I,6) 

0069400 

C. . . . 

.LAYER STIFFNESS MATRICES C 

0069500 


CALL SKINV(ST,6,SF,6,VC6,6,-1,SIMUL,6) 

0069600 


URITE(IN6,112) K 

0069700 

112 

FORMAT (/IX, 'STIFFNESS MATRIX C FOR ',12,' MATERIAL') 

0069800 


DO 108 1=1,6 

0069900 

108 

HRITE(IN6,123)SF(I,1),SF(I,2).SF(I,3),SF(I,4),SF(I,5>,SF(I,6> 

0070000 

C • • • • 

.PLANESTRESS STIFFNESS MATRIX 

0070100 


CALL PZERO (T2, 36) 

0070200 


CALL PZERO (Tl, 36) 

0070300 


T1(1,1)=ST(1,1) 

0070400 


T1(1,2)=ST(1,2) 

0070500 


T1(1,3)=ST(1,6) 

0070600 


T1(2,1)=ST(2,1) 

0070700 


T1(2,2)=ST(2,2) 

0070800 


T1(2,3)=ST(2,6) 

0070900 


T1(3,1)=ST(6,1) 

0071000 


Tl (3,2)=ST(6,2) 

0071100 


T1(3,3)=ST(6,6) 

0071200 


Tl (4,4)=1. 
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0071300 

00714100 

0071500 

0071600 

0071700 

0071800 

0071900 

0072000 

0072100 

0072200 

0072300 

0072400 

0072500 

0072600 

0072700 

0072800 

0072900 

0073000 

0073100- 

0073200 

0073300 

0073400 

0073500 

0073600 

0073700 

0073800 

0073900 

0074000 

0074100 

0074200 

0074300 

0074400 

0074500 

0074600 

0074700 

0074800 

0074900 

0075000 

0075100 

0075200 

0075300 

0075400 

0075500 

0075600 

0075700 

0075800 

0075900 

0076000 

0076100 

0076200 

0076300 

0076400 

0076500 

0076600 


T1(5,5>=1. 

Tl(6,6)=l. 

CALL SKINV(T1,6,T2,6,VC6,6,-1,SIMUL,6) 

CALL PZER0(T1,36) 

T1(1,1)=-ST(1,3) 

T1(2,1)=-ST(2,3) 

T1(3,1)=-ST(6,3) 

CALL BID0T(T2,6,T1,6,T1I,6) 

WRITE (IN6, 115) 

115 FORMAT ( /IX , 'PLANESTRESS STIFFNESS MATRIX CBB r CB OF M2,' LAYER') 

DO 113 1=1,3 

113 HRITE(IN6,123) T2(I,1),T2(I,2),T2(I,3),T1I(I,1) 

C MATERIAL PROPERTIES FOR PLY LEVEL ANALYSIS 

S31(K)=ST(3,1> 

S32(K)=ST(3,2) 

S33(K)=ST(3,3) 

S44(K)=ST(4,4) 

S45(K)=ST(4,5) 

S55(K)=ST(5,5) 

CB1 (K)=T1I (1,1 ) 

CB2(K>=T1I(2,1) 

CB6(K)=T1I(3,3) 

CB11(K)=T2(1,1) 

CB12(K)=T2(1,2) 

CB16(K)=T2(1,3) 

CB22(*')=T2(2,2> 

CB26(K)=T2(2,3) 

CB66(K)=T2(3,3) 

SNC1(K)=CB12(K)/CB22(K)*(S31(K>*CB12(K>+S32(K>*CB22(K)>+CB12(K)*S4- 

•4(K)-CB1(K) 

SNC2(K)=CB22(K)/CB22<K)*(S31(K)»CB12(K)+S32(K)*CB22(K))+CB22(K)*S4- 

•4(K)-CB2(K) 

KPl(K)=CBl(K>/2 
KP2(K)=CB2(K)/2 
SHM2 (K ) =S44 (K>* ( - .25 > 

SHQ2(K)=S44(K)*.75 

KQ1(K)=(.25*CB12(K)*(S31(K>*CB12(K>/CB22(K) +S32(K> )-CB12(K)»S44(K- 
*>+CBl(k>>*.6 

KQ2(K)=(.25*CB22(K)*(S31(K)*CB12(K)/CB22(K) +S32(K) >-CB22(K)*S44(K- 
*>+CB2(K>>*.6 

KM1(K)=(1.5*CB12(K)*(S31(K)*CB12(K)/CB22(K) +S32(K))-CB12(K)*S44(K- 
•)+CBl(K))*.l 

KM2(K>=(1.5*CB22(K)*(S31(K)*CB12(K)/CB22(K> +S32(K> >-CB22(K)*S44(K- 
•)+CB2(K))*.l 

KN1(K)=(CB12(K)*(S31(K>*CB12(K)/CB22(K) +S32(K) >+CB12(K)»S44(K>+2*- 
•CB1(K))/12 

101 KN2(K>=(CB22(K>*(S31(K>*CB12(K)/CB22(K> +S32(K> >+CB22(K)*S44(K>+2*- 

•CB2(K))/12 
WRITE (IN6, 125) 

125 FORMATC/1X, r CBll,CB12,CB16,CB22,CB26,CB66,CBl,CB2,CB6 r ) 

123 FORHAT(1X,8E10.4) 

DO 126 K=1,N 

126 WRITE ( IN6 , 123 )CB11(K) ,CB12(K) ,CB16(K)»CB22(K) »CB26(K> ,CB66(K) , 
*CB1(K),CB2(K),CB6(K> 
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0076700 

0076800 

0076900 

0077000 

0077100 

0077200 

0077300 

0077300 

0077500 

0077600 

0077700 

0077800 

0077900 

0078000 

0078100 

0078200 

0078300 

0078300 

0078500 

0078600 

0078700 

0078800 

0078900 

0079000 

0079100 

0079200 

0079300 

0079300 

0079500 

0079600 

0079700 

0079800 

0079900 

0080000 

0080100 

0080200 

0080300 

0080300 

0080500 

0080600 

0080700 

0080800 

0080900 

0081000 

0081100 

0081200 

0081300 

0081300 

0081500 

0081600 

0081700 

0081800 

0081900 

0082000 


URITE C IN6 , 129 ) 

1 29 FORMAT ( / IX , ' S31 , S32 , S33 f S33 , S35 , S55 , SNC1 , SNC2 , SHM2 , SHQ2 r ) 

DO 130 K=1,N 

130 URITE (IN6, 123 )S31 (K) ,S32(K) ,S33(K),S33(K),S35(K),S55(K),SNC1 (K), 
*SNC2(K),SHM2(K),SHQ2(K) 

URITE(IN6,119> 

119 F0RMAT(/1X,'KN1,KN2,KP1,KP2,KM1,KM2,KQ1,KQ2'> 

DO 122 K=1,N 

122 URITE (IN6, 123) KN1(K),KN2(K),KP1(K),KP2(K),KM1(K>,KM2(K),KQ1(K), 
•KQ2(K> 

RETURN 

END 

FUNCTION UVD(N,I,Z,J,RCM,IPR) 

IMPLICIT REAL <A-H,0-Z) 

REAL RCM 
DIMENSION RCH(l) 

M=«J*IPR 

UVD= UVD1 (N,I,Z, RCM ( 1 ) , RCM ( 1+M ) , RCM ( 1+2*M> , RCM ( 1+3*M) , RCM( 1+3*M ) 
• , RCM ( 1+5*M ) , RCM ( 1+6*M ) , RCM ( 1+7*M ) , RCM ( 1+8*M> , RCM( 1+9»M ) 
•,RCM(1+10*M>,RCM(1+11*M),RCM(1+12*M>,RCM(1+13*M),RCM(H-13*M) 
•,RCM(1+15*M) ,RCM(1+16*M) ,RCM(1+17*M) ,RCM(1+18*M),RCH(1+19*M) 

* , RCM ( 1+20*M ) , RCM ( 1 +21 *M ) , RCM ( 1+22*M ) , RCM ( 1+23*M > , RCM O +23*H ) 

* , RCM ( 1+25*M) , RCM ( 1+26*M) , RCM ( 1+27*M> , RCM ( 1+28*M) , RCM( 1+29*M > 
•,RCM(l+30*M),RCM(l+31*M),RCM(l+32*M),RCM<l+33*M> r RCMU+33*M) 

• , RCM ( 1+35*M ) , RCM ( 1+36*M ) , RCM ( 1+37 *M ) , RCM ( 1+39*M ) r RCM ( 2+39*M ) ) 
RETURN 
END 

FUNCTION UVD1(N,I,Z, Ell,E22,E33,NU12,NU13,NU23,G33,G55,G66 r S31, 
*S32 , S33 r S33 , S35 , S55 , CB1 , CB2 , CB6 , CB1 1 , CB12 , CB16 r CB22 , CB26 , CB66 , SNC1 
* , SNC2 , SHM1 r SHM2 , SHQ1 , SHQ2 , KQ1 , KQ2 , KM1 , KM2 , KP1 , KP2, KN1 , KN2, EPSLON , 
•CURVTR) 

IMPLICIT REAL <A-H r O-Z> 

REAL KMl,KM2,KNl,KN2,KQl,KQ2,KPl f KP2,NU12,NU13 r NU23 
DIMENSION Ell(l),E22(l) r E33(l),NU12(l),NU13(l),NU23(l),G33(l), 
*G55(l>,G66(l) r S31(l),S32(l),S33(l),S33(l),S35(l),S55(l),CBl(l), 
•CB2(l),CB6(l),CBll(l),CB12(l),CB16(l),CB22(l),CB26(l) f CB6r>(l), 

*SNC 1(1), SNC2 ( 1 ) , SHM1 ( 1 ) , SHM2 ( 1 ) » SHQ1 ( 1 ) , SHQ2 ( 1 ) , KQ1 ( 1 ) , KQ2 ( 1 ) , 

•KM1 (1 ),KM2(1 ) ,KP1 (1 ) ,KP2(1 ) f KN1 (1 ) ,KN2(1 ) 

GOTO (1,2, 3, 3 ,5, 6, 7, 8, 9, 10, 11, 12, 13, 13, 15, 16, 17, 18), N 

1 UVD1=1 
RETURN 

2 UVD1 =- (S31 ( I )*CB12 ( I )+S32 ( I )*CB22 ( I ) )/2*Z»Z 

RETURN 

3 UVD1 = ( S31 ( I ) »CB1 1(1) +S32 ( I ) *CB12 ( I ) )* (Z*EPSL0N-Z*Z/2*CURVTR ) 

C HERE Z SHOULD HAVE A DIMENSION OF A LENGTH AND U,XX SHOULD 

C BE GIVEN AS IS UITHOUT ANY NONDIMENSION ALIZATION 

RETURN 

3 UVD1 =<S31(I)*CB12(I)+S32(I)*CB22(I))*Z 

RETURN 

5 UVD1 =(S31(I)*SNC1(I)+S32(I)*SNC2(I))/12«(Z*Z*Z-Z) 

•-S33 ( I )/12* (Z**3-3*Z )+ (S31 ( I >*KN1 ( I )+S32 ( I )*KN2 ( I ) )*Z 

RETURN 

6 UVD1 =(S31(I)*SNCl(I)+S32(I)*SNC2(I))/3*(Z**3/3-3./10*Z*Z) 

•-S33 ( I )/3* (Z*»3/3-Z*Z/2 )+ (S31 ( I >*KM1 ( I )+S32( I )*KM2(I ) )/2*Z»Z 
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0082100 


0082200 

7 

0082300 


00824100 

8 

0082500 


0082600 


0082700 

9 

0082800 


0082900 

10 

0083000 


0083100 


0083200 

11 

0083300 


0083000 

12 

0083500 


0083600 

13 

0083700 


0083800 

14 

0083900 


0080000 


0080100 


0080200 

15 

0080300 


0080000 

16 

0080500 


0080600 


0080700 


0080800 


0080900 


0085000 

17 

0085100 


0085200 

18 

0085300 


0085000 


0085500 


0085600 


0085700 


0085800 


0085900 


0086000 


0086100 


0086200 


0086300 


0086400 


0086500 


0086600 


0086700 


0086800 


0086900 


0087000 

80 

0087100 


0087200 

8 

0087300 


0087400 

2 


RETURN 

HVD1 = (S31 ( I )*KP1 ( I >+S32 ( I >*KP2 ( I > ) *Z+S33 ( I ) *Z/2 
RETURN 

UVD1 = (S31 ( I )*KQ1 ( I )+S32 ( I >*KQ2 (I ) ) /2*Z*Z-S?3 < I >/4l» (Z**4/4-Z*Z/2 ) 
*+ ( S31 ( I >*SNC1 ( I ) +S32 ( I > *SNC2 ( I ) ) /4* ( Z**4/4-Z*Z*3/l 0 ) +S33 (I)/4*Z*Z 
RETURN 
UVD1 =-Z 

RETURN 

UVD1 = (S31 ( I > *CB12 ( I ) +S32 ( I ) *CB22 CD) /6*Z**3+S44 ( I ) *CB22 (I)/6*( 

•Z**3-3*Z> 

RETURN 
HVD1 =1 

RETURN 

WVD1 =-S31(I)*CB12(I)/2*Z*Z 

RETURN 

HVD1 =<S44(I>+S32(I>)/4*Z*Z 

RETURN 

WVD1 =-S31(I)*KNl(I)/2*Z*Z 

*- ( S31 ( I )*SNC1 ( I > +S32 ( I ) *SNC2 (I))/12*( Z**4/4 - 2 * 2 / 2 ) 

•+S33CI )/12* (Z**4/4-3. /2*Z*Z )+0*S44 C I )/4*SNC2 ( I )* ( -Z**4/12+Z**2/6 ) 
RETURN 
HVD1 
RETURN 
HVD1 

• 

• 

RETURN 
UVD1 
RETURN 
WVD1 


RETURN 
END 

SUBROUTINE CDET ( NSYM , COEF, X,F, NSDFG , CVC30 ,NPIVOT,NK10, 

•N , NRT ,NCT,CRT, NRGN ,NT,NG2, SKO , SKI ,SK2,SK3,SK4, A, WA,NCN1, N25L ) 
IMPLICIT REAL <A-H,0-Z> 

COMPLEX A(N25L, 1 ) , CVC30 ( 1 > ,X, F, X2347, DCMPLX 

DIMENSION WA(1 > ,COEF(l > ,SK0 (N25L, 1 > ,SK1 (N25L, 1 ) ,SK2(N25L, 1 ) , 
•SK3CN25L, 1 ) ,SK4 (N25L, 1 > ,CRT (NCN1 , 1 ) 

IF (NSDFG. NE.O) GOTO 8 
IF(NKIO.EQ.O) GOTO 8 
F=<0.D 00,0.D 00) 

DO 80 1=1, NT 
F=F+COEF ( I )*X** (2*1-2 ) 

RETURN 
DO 2 1=1, N 
DO 2 J=1,N 

A(I, J)=SK0(I, J)+SK1 (I, J)*X+SK2(I, J)*X*X +SK3(I, J)*X*X*X 


=S44(I)/2*Z 

=- (S31 ( I )*KM1 ( I)+S32 ( I )*KM2 ( I ) >/6*Z**3 
-S44(I)*KM2(I)/6*(Z**3-3*Z> 

-(S31 (I )*SNC1 (I )+S32(I )*SNC2(I ) )/4* (Z**5/20-Z**3/10 ) 
+S33(I)/4*(Z**5/20-Z**3/6> 

+0*S44 ( I ) /4»SNC2 ( I )* ( -Z**5/20+Z**3/l 0-Z/20 ) 

=-S31 ( I )*KP1 (I )*Z*Z/2-S33 ( I )*Z*Z/4 

=-(S31(I)*KQl(I)+S32(I)*KQ2(I))/6*Z**3 
-S44 ( I )*KQ2 ( I )/6* (Z**3-3*Z ) 

- <S31 ( I )*SNC1 (I)+S32 ( I )*SNC2( I ) >/4* (Z**5/20-Z**3/10 ) 
-S33 ( I ) *Z**3/12 
♦S33 ( I ) * ( Z**5/20 - Z** 3/6 > /4 
♦0*S44 ( I ) /4*SNC2 ( I )* ( -Z**5/20+Z**3/l 0-Z/20 ) 
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0087500 

0087600 

0087700 

0087800 

0087900 

0088000 

0088100 

0088200 

0088300 

0088400 

0088500 

0088600 

0088700 

0088800 

0088910 

0089000 

0089100 

0089200 

0089300 

0089400 

0089500 

0089600 

0089700 

0089800 

0089900 

0090000 

0090100 

0090200 

0090300 

0090400 

0090500 

0090600 

0090700 

0090800 

0090900 

0091000 

0091100 

0091200 

0091300 

0091400 

0091500 

0091600 

0091700 

0091800 

0091900 

0092000 

0092100 

0092200 

0092300 

0092400 

0092500 

0092600 

0092700 

0092800 


*+SK4(I,J)»X*X*X+X 
IF (NSDFG.NE.O) GOTO 3 

CALL LEQT1 C ( A , N , N25L , CV C3 0 , 1 , N25L , 0 , HA , I ER ) 

F=(1.D 00,0. D 00) 

DO 1 1=1, N 
IPVT=WA(I> 

IF(IPVT.NE.I) F=-F 
1 F=F*A(I,I) 

IF(NPIVOT.EQ. 11000) RETURN 
Y=1.D 00 
NRTNCT =NRT +NCT 
IF (NRTNCT.LT. 2) GOTO 232 
11=1 

215 IF(CRT(NG2+I1,2).NE.0. ) GOTO 231 

Y=Y<* (X**2-CRT(NG2+I1 , 1 )**2 ) 

11=11+2 
GOTO 214 

231 Y=Y* ( (X-CRT (NG2+I1 , 1 ) )**2+CRT (NG2+I1 , 2 )**2 ) 

Y=Y* ( <X+CRT(NG2+I1 , 1 ) )**2+CRT (NG2+I1 ,2 )**2) 

11=11+4 

214 IFdl.GT. NRTNCT) GOTO 232 
GOTO 215 

232 F=F/Y 
IF(NSYM.EQ.l) F=F/X**5 
IF (NSYM.EQ.2) F=F/X**3 
IF'N’SYM.EQ.3> F=F/X**6 
RET'oftN 

3 CVC30(NPIVOT)=(1.D 00,0.D 00) 

A(NPIV0T,NPIVOT)=(l.D 00,0.D 00) 

DO 4 1=1, N 

IF ( I . NE . NPI VOT ) CVC30 ( I ) =- A ( I , NPI VOT ) 

IF (I. HZ. NPI VOT) A(NPIVOT, I ) = (0.D 00,0.D 00) 

4 IF ( I. NE. NPI VOT) A(I,NPIVOT)=(0.D 00,0.D 00) 

CALL LEQT1C ( A, N , N25L , CVC30 , 1 , N25L , 0 , MA , IER ) 

F=(1.D 00, O.D 00) 

DO 71 1=1, N 
IPVT=HA(I) 

IF(IPVT.NE.I) F=-F 
71 F=F*A(I,I> 

RETURN 

END 

SUBROUTINE FDET (NSYM, COEF, X , F ,NSDFG, VC30 , NPIVOT, NK10 ,NCN1 , N, NRT , 
•NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A , N25L ) 

IMPLICIT REAL (A-H,0-Z) 

DIMENSION CRT(NCN1 , 1 ) ,SK0 (N25L, 1 ) ,SK1 (N25L, 1 ) ,SK2(N25L, 1 ) , 

•SK3 (N25L, 1 ) »SK4 (N25L , 1 ) , A(N25L , 1 ) , COEF ( 1 ) , VC30 ( 1 ) 

IF(NSDFG.NE.O) GOTO 8 
IF(NKIO.EQ.O) GOTO 8 
F=0.D 00 
DO 80 1=1, NT 

80 F=F+C0EF(I)»X**(2*I-2) 

RETURN 

8 DO 2 1=1, N 

DO 2 J=1,N 
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0092900 2 

0093000 

0093100 

0093200 

0093300 

0093400 

0093500 

0093600 

0093700 

0093800 215 

0093900 

0094000 

0094100 

0094200 231 

0094300 

0094400 

0094500 214 

0094600 

0094700 '232 

0094800 

0094900 

0095000 

0095100 

0095200 3 

0095300 

0095400 

0095500 

0095600 

0095700 4 

0095800 

0095900 

0096000 

0096100 

0096200 

0096300 

0096400 

0096500 

0096600 

0096700 

0096800 

0096900 

0097000 

0097100 

0097200 

0097300 

0097400 

0097500 

0097600 

0097700 215 

0097800 

0097900 

0098000 

0098100 231 

0098200 


A(I,J)=SK0(I,J)+SK1(I,J)*X+SK2(I,J)*X*X+SK3(I,J)*X*X*X 
•+SK4(I,J)*X*X*X*X 
IF (NSDFG.NE.O) GOTO 3 
CALL SKINV(A,N25L,A,N25L,VC30,1,-1,F,N) 

IF(NPIVOT.EQ. 11000) RETURN 
Y=1.D 00 
NRTNCT =NRT+NCT 
IF (NRTNCT.LT. 2) GOTO 232 
11=1 

IF(CRT(NG2+I1,2).NE.0. ) GOTO 231 
Y=Y* (X**2-CRT < NG2+I1 , 1 )**2 ) 

11=11+2 
GOTO 214 

Y=Y* ( ( X-CRT < NG2+I1 , 1 ) )**2+CRT ( NG2+1 1 , 2 )**2 ) 

Y=Y* ( (X+CRT (NG2+I1 , 1 ) )**2+CRT (NG2+I1 ,2>**2 ) 

11=11+4 

IF ( II. GT. NRTNCT) GOTO 232 

GOTO 215 

F=F/Y 

IF (NSYM.EQ. 1 ) F=F/X**5 
IF (NSYM.EQ. 2) F=F/X*»3 
IF (NSYM.EQ. 3) F=F/X**6 
RETURN 

VC30(NPIVOT)=1.D 00 
A(NPIVOT,NPIVOT)=l.D 00 
DO 4 1=1, N 

IF ( I . NE . NPI VOT ) VC30 ( I ) =- A ( I , NPI VOT ) 

IF(I.NE.NPIVOT) A(NPIVOT,I )=0.D 00 
IF ( I. NE. NPI VOT) A(I,NPIVOT)=O.D 00 
CALL SKINV ( A, N25L, A, N25L, VC30 , 1 , 0 , F, N) 

RETURN 

END 

SUBROUTINE COEFF1 (LJ , COEF , N , NCT , NRT , NRGN, NT, NG2, SK, SKO , 

•SKI , SK2 , SK3 , SK4 , CRT , NK10 , NCN1 , A, N25L ) 

IMPLICIT REAL (A-H,0-Z) 

DIMENSION SK(N25L, 1 ) , C0EF(1 ) ,SK0 (N25L, 1 > ,SK1 (N25L, 1 ) ,SK2(N25L, 1 ) , 
•SK3(N25L, 1 ) ,SK4 (N25L, 1 ) ,CRT (NCN1 , 1 ) , A(N25L, 1 ) 

DO 6 1=1, NT 
Y=DFLOAT(I) 

X=1.D 00 

IF (LJ.EQ.2) X=X*Y*Y*Y 
IF(LJ.EQ.l) X=X*Y*Y*Y*Y*Y 
IF(LJ.EQ.3) X=X*Y*Y*Y*Y*Y*Y 

CALL FDET(LJ, COEF, Y, FI, 0, COEF, 11000, NK10,NCN1, 

•N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A, N25L ) 

NRTNCT =NRT +NCT 

IF (NRTNCT.LT. 2) GOTO 232 

11=1 

IF (CRT(NG2+I1,2) .NE.O . ) GOTO 231 
X=X« ( Y**2-CRT(NG2+I1, 1 )**2 ) 

11=11+2 
GOTO 214 

X=X* ( ( Y-CRT ( NG2+I1 , 1 ) ) »»2+CRT ( NG2+I 1 , 2 )«»2 ) 

X=X*(( Y+CRT ( NG2+I1 , 1 ) ) **2+CRT ( NG2+I 1 ,2)**2) 
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0098300 

0098400 

0098500 

0098600 

0098700 

0098800 

0098900 

0099000 

0099100 

0099200, 

0099300 

0099400 

0099500 

0099600 

0099700 

0099800 

0099900 

0100000 

0100100 

0100200 

0100300 

0100400 

0100500 

0100600 

0100700 

0100800 

0100900 

0101000 

0101100 

0101200 

0101300 

0101400 

0101500 

0101600 

0101700 

0101800 

0101900 

0102000 

0102100 

0102200 

0102300 

0102400 

0102500 

0102600 

0102700 

0102800 

0102900 

0103000 

0103100 

0103200 

0103300 

0103400 

0103500 

0103600 


11=11+4 

214 IF(Il.GT.NRTNCT) GOTO 232 
GOTO 215 

232 COEF(I)=FI/X 
DO 6 J=1,NT 

6 SK(I,J)= Y**(2*J-2) 

CALL SKINV(SK,N25L,SK,N25L» COEF , 1 , 0 , SIMUB , NT ) 

RETURN 

END 

SUBROUTINE CROOT ( NSYtt , N , NT , NG2 , NRGN , NK1 0 , NRT , NCT , IN7 , IN8 , 

•CRT , CMODE, SKO , SKI , SK2 , SK3 , SK4 , COEF , CVC31 , CVC32 , A, MA, NCN1 , NCN2 , N25L 

• ) 

IMPLICIT REAL (A-H,0-Z) 

DIMENSION CRT(NCN1 , 1 ) ,CM0DE(NCN1 ,NCN2, 1 ) ,SK0 (N25L, 1 ) ,SK1 (N25L, 1 ) , 
•SK2(N25L,1),SK3( N25L , 1 ) , SK4 ( N25L , 1 > , COEF ( 1 ) , MA ( 1 ) 

COMPLEX XF12,FF1,FF2,XI,XF,FI,CVC31(1),CVC32(1),XRL,XIM, 
•A(N25L, 1 ) , DCHPLX 

C LIMIT MAX. NO. OF ITERATIONS TO 100, AND ROOTS TO 100 

XRL=(1.,0.) 

XIM=(0.,-1.) 

DO 4 1=1,100 
4011 HRITE(IN8,5> 

5 F0RMAT(1X, r 2F10.6 XI.. INPUT INITIAL COMPLEX ROOT VALUE') 

READ(IN7,6) XI 

IF (CABS (XI ) .LE.l.D-20 ) GOTO 14 

6 FORMAT (10F10. 6) 

NITER=0 

9 CALL CDET(NSYM,COEF,XI,FI,0,CVC31,1»NK10, 

•N,NRT,NCT,CRT,NRGN,NT,NG2,SK0,SK1,SK2,SK3,SK4, A,HA,NCN1,N25L) 
RFI=FI*XRL 
RGI=FI*XIM 
UI=XI*XRL 
VI=XI*XIM 

XI=CMPLX( 1.00000000001*UI,VI) 

CALL CDET(NSYM,COEF,XI,FI,0,CVC31,1,NK10, 

•N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A , WA , NCN1 , N25L ) 
RFU=FI*XRL 
RGU=FI<*XIM 

RFU=(RFU-RFI)/.00000000001/UI 

RGU=(RGU-RGI)/.00000000001/UI 

XI=CMPLX(UI,1.00000000001*VI) 

CALL CDET(NSYM,COEF,XI,FI,0,CVC31,1,NK10, 

•N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A , WA , NCN1 , N25L ) 
RFV=FI*XRL 
RGV=FI*XIM 

RFV=(RFV-RFI)/.00000000001/VI 
RGV=(RGV-RGI )/.00000000001/VI 
UF=UI- (RFI»RGV-RFV*RGI ) / (RFU*RGV-RFV*RGU> 

VF=VI- (RFU*RGI-RFI*RGU )/ (RFU«RGV-RFV*RGU ) 

NITER=NITER+1 
URITECIN8 , 12 )NITER, XI , FI 

IF(CABS(XI/CMPLX(UF,VF)-1. ).LE. .000000001) GOTO 7 

IF (NITER. GT. 100) GOTO 4011 

XI=CMPLX(UF,VF) 
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0103700 GOTO 9 

0103800 7 WRITE ( IN8 , 12 ) I , XI , FI 

0103900 WRITE(IN8,76> 

0104000 READ(IN7,17) J 

0104100 17 FORMAT (2012) 

0104200 IF(J.NE.O) GOTO 77 

0104300 GOTO 4011 

0104400 12 FORMAT (IX, 13, IX, 4D15. 8, r I, XI, FI') 

0104500 77 DO 7712 1111=1,4 

0104600 CRT (NG2+NRT+NCT+IIII,1 )=UF 

0104700 7712 CRT(NG2+NRT+NCT+IIII,2)=VF 

0104800 NCT=NCT+4 

0104900 76 FORMAT (IX, '12... INPUT... J=1 TO STORE ROOT') 

0105000 4 CONTINUE 

010 SlOO 14 WRITE(IN8,21) 

0105200 21 FORMAT (IX, '12... INPUT J=0 TO FIND COMPLEX EIGEN MODES') 

0105300 READ(IN7, 17) J 

0105400 IF(J.NE.O) GOTO 66 

0105500“ 111=0 

0105600 DO 82 I=1,NCT,4 

0105700 111=111+1 

0105800 X1=CRT (NG2+I , 1 ) 

0105900 X2=CRT(NG2+I,2) 

0106000 DO 821 Jl=l,4,2 

0106100 IF(Jl.NE.l) GOTO 832 

0106200 XI =CMPLX ( XI, X2) 

0106300 CALL CDET ( NSYM, COEF , XI , FI , 1 , CVC31 , 1 , NK1 0 , 

0106400 *N,NRT,NCT,CRT,NRGN,NT,NG2,SK0,SK1,SK2,SK3,SK4,A,WA,NCN1,N25L) 

0106500 861 FORMAT (/3 (IX, 14), IX, 2D15. 8, IX, 'COMPLEX EIGEN VALUE AND MODE') 

0106600 832 IF(Jl.NE.3) GOTO 821 

0106700 XI=CMPLX(-X1, X2) 

0106800 CALL CDET(NSYM,COEF,XI,FI,1,CVC32,1,NK10, 

0106900 •N,NRT,NCT,CRT,NRGN,NT,NG2,SK0,SK1,SK2,SK3,SK4,A,WA,NCN1,N25L) 

0107000 87 FORMAT ( 10 ( IX, D10.4 ) ) 

0107100 821 CONTINUE 
0107200 NCT12=NCT/2 

0107300 IF(NRGN.EQ.O) WRITE (IN8, 861 ) NG2,NCT,I,XI 

0107400 If (NRGN.EQ.l) WRITE(IN8,861) NG2,NCT12,III,XI 

0107500 IF(NRGN.EQ.l) CRT(NG2+2*III-1,1)=X1 

0107600 IF(NRGN.EQ.l) CRT(NG2+2*III-1,2)=X2 

0107700 IF(NRGN.EQ.l) CRT(NG2+2*III-0,1)=X1 

0107800 IF (NRGN.EQ.l) CRT(NG2+2*III-0,2)=X2 

0107900 DO 7634 J=1,N 

0108000 7634 WRITE(IN8,87) CVC31(J),CVC32( J) 

0108100 DO 823 J=1,N 

0108200 RCVC3 1 =CVC31 ( J ) *XRL 

0108300 ACVC31=CVC31(J)*XIM 

0108400 RCVC32=CVC32(J)*XRL 

0108500 ACVC32=CVC32(J)*XIM 

0108600 IF (NRGN.EQ.l) GOTO 76341 

0108700 CMODE ( NG2+I+0 , J , 1 ) = ( ACVC31+ACVC32)/2. 

0108800 CMODE(NG2+I+0,J,2)=( RCVC31-RCVC32)/2. 

0108900 CMODE ( NG2+I+0 , J , 3 ) = ( ACVC31~ACVC32)/2. 

0109000 CMODE (NG2+I+0 , 0 , 4 ) = ( RCVC31+RCVC32)/2. 
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0109100 CHODE ( NG2+I+1 , J , 1 ) = ( ACVC31-ACVC32)/2. 

0109200 CH0DE(NG2+I+1,J,2)=( RCVC31+RCVC32)/2. 

0109300 CHODE ( NG2+I+1 , J , 3 ) = ( ACVC31+ACVC32)/2. 

0109400 CHODE ( NG2+I+1 , J , 4 ) = ( RCVC31-RCVC32>/2. - 

01C9500 CH0DE(NG2+I+2,0,1)=( RCVC31+RCVC32)/2. 

0109600 CHODE (NG2+I+2 , J , 2 ) = ( - ACVC31+ACVC32 ) /2 . 

0109700 CHODE (NG2+I+2, J,3)=( RCVC31-RCVC32)/2. 

0109800 CHODE ( NG2+I+2, J , 4 ) = ( - ACVC31 - ACVC32 ) /2 . 

0109900 CHODE (NG2+I+3, J,l)=( RCVC31-RCVC32)/2. 

0110000 CH0DE(NG2+I+3,J,2)=(-ACVC31-ACVC32)/2. 

0110100 CH0DE(NG2+I+3,<J,3) = ( RCVC31+RCVC32)/2. 

0110200 CHODE ( NG2+I+3 , J , 4 > = < -ACVC31+ACVC32 )/2. 

0110300 GOTO 823 

0110400 76341 CONTINUE 

0110500 CH0DE(NG2+2*III-1,J,1>= RCVC31 

0110600 CHODE ( NG2+2* III-1,J,2) =- ACVC31 

0110700 CHODE(NG2+2*III-1,J,3)=O.D 00 

0110800 CHODE(NG2+2*III-1,J,4)=O.D 00 

0110900' CHODE(NG2+2*III-0, J,l> = ACVC31 

0111000 CHODE(NG2+2*III-0,J,2)= RCVC31 

0111100 CHODE (NG2+2*III-0,J,3)=0.D 00 

0111200 CHODE (NG2+2*III-0, J r 4)=0.D 00 

0111300 823 CONTINUE 

0111400 82 CONTINUE 

0111500 IF(NRGN.EQ.l) NCT=NCT/2 

0111600 86 RETURN 

0111700 END 

0111800 SUBROUTINE RROOT ( NSYH , N , NT f NG2 , NRGN , NK1 0 , NRT , NCT , IN7 , 1 N8 , 

0111900 *CRT , CHODE , SKO , SKI , SK2 , SK3 , SK4 , COEF , VC31 , VC32 , A, NCN1 - NCN2 , N25L ) 

0112000 IHPLICIT REAL (A-H,0-Z) 

0112100 DIHENSION CRTCNCNl > l),CH0DE(NCNl,NCN2,l),SK0(N25L > l) f SKl(N25L r l), - 

0112200 *SK2 ( N25L , 1 > , SK3 ( N25L , 1 ) , SK4 ( N25L , 1 ) , COEF (1)»VC31(1), VC32 < 1 ) 

0112300 *, A(N25L r l > 

0112400 C6390 F0RJ1AT(1X,I2,14F9.4) 

0112500 C NQQ=N 

0112600 C DO 6391 1=1, NQQ 

0112700 C6391 HRITE(6,6390 ) I, (SKO (I, J) , J=1,NQQ> 

0112800 C DO 6392 1=1, NQQ 

0112900 C6392 MRITE(6,6390 ) I, <SK1(I, J),J=1,NQQ> 

0113000 C DO 6393 1=1, NQQ 

0113100 C6393 WRITE(6,6390> I, (SK2(I, J> , J=1,NQQ> 

0113200 C DO 6394 1=1, NQQ 

0113300 C6394 WRITE(6,6390 ) I, (SK3CI, J),J=1,NQQ) 

0113400 C DO 6395 1=1, NQQ 

0113500 C6395 WRITE(6,6390> I, (SK4(I, J) , J=1,NQQ) 

0113600 C LIHIT HAX. NO. OF ROOTS TO 100, AND ITERATIONS TO 100 

0113700 DO 4 1=1,100 

0113800 4011 HRITE(IN8,5) 

0113900 5 FORMAT C IX, r 10F10.6 XI.. INPUT INITIAL ROOT VALUE') 

0114000 READ(IN7,6> XI 

0114100 IF(DABS(XI).LE.l.D-20) GOTO 14 

0114200 6 FORMAT (10F1 0.6) 

0114300 NITER=0 

0114400 9 CALL FDET(NSYH,COEF,XI,FI,0,VC31,1,NK10,NCN1, 
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0114500 

0114600 

0114700 

0114800 

0114900 

0115000 

0115100 

0115200 

0115300 

0115400 

0115500 

0115600 

0115700 

0115800 

0115900 

0116000 

0116100 

0116200 

0116300 

0116400 

0116500 

0116600 

0116700 

0116800 

0116900 

0117000 

0117100 

0117200 

0117300 

0117400 

0117500 

0117600 

0117700 

0117800 

0117900 

0118000 

0118100 

0118200 

0118300 

0118400 

0118500 

0118600 

0118700 

0118800 

0118900 

0119000 

0119100 

0119200 

0119300 

0119400 

0119500 

0119600 

0119700 

0119800 


*N , NRT , NCI , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A , N25L ) 
XF1=(1.0+0.00000000001)*XI 

CALL FDET (NSYM, COEF, XF1 , FF1 , 0 , VC31 , 1 ,NK10,NCN1 , 

•N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A , N25L ) 
XF2=(1. 0-0. 00000000001 )*XI 

CALL FDET ( NSYM , COEF ,XF2,FF2,0,VC31,1,NK10,NCN1, 

•N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A, N2SL ) 
DERXI = ( FF1 -FF2 ) / ( XF1 -XF2 ) 

XF=XI-FI/DERXI 

NITER=NITER+1 

HRITE ( IN8 , 12 )NITER r XI , FI 

IF(DABS(XI/XF-1.).LE.. 000000001) GOTO 7 

IF (NITER. GT . 100 ) GOTO 4011 

XI=XF 

GOTO 9 

7 WRITE ( IN8 , 12 ) I , XI , FI 

12 FORMAT(lX,I3,lX,2(D15.8,lX) r ’ I, XI, FI') 

WRITE (IN8, 76) 

READ(IN7,17) J 
17 FORMAT (2012) 

IF(J.NE.O) GOTO 77 
GOTO 4011 

77 DO 7712 1111=1,2 

CRT(NG2+NRT+NCT+IIII,2)=0.D 00 
7712 CRT(NG2+NRT+NCT+IIII,1)=XI 
NRT=NRT+2 

76 FORMAT (IX, '12... INPUT... J=1 TO STORE ROOT') 

4 CONTINUE 

14 WRITE (INS, 21) 

21 FORMAT (IX, '12... INPUT J=0 TO FIND REAL EIGEN MODES') 
READ(IN7,17) J 
IF(J.NE.O) GOTO 66 
111=0 

DO 82 1=1, NRT, 2 
111=111+1 
XI=-CRT(NG2+I,1) 

CALL FDET (NSYM, COEF, XI,FI,1,VC32,1,NK10,NCN1, 

•N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A, N25L ) 
XI=CRT (NG2+I , 1 ) 

CALL FDET ( NSYM, COEF, XI , FI , 1 , VC31 , 1 , NK10 , NCN1 , 

*N , NRT , NCT , CRT , NRGN , NT , NG2 , SKO , SKI , SK2 , SK3 , SK4 , A, N25L ) 
DO 83 11=1, N 
IF(NRGN.EQ.l) GOTO 24850 
CM0DE(NG2+I ,Il,l)=(VC31(Il)+VC32(Il))/2 
CM0DE(NG2+I ,Il,2)=(VC31(Il)-VC32(Il))/2 
CM0DE(NG2+I+l,Il,l)=(VC31(Il)-VC32(Il))/2 
CMODE ( NG2+I+1 , I 1 , 2 ) = ( VC31 ( I 1 ) + VC32 ( I 1 ) ) /2 
GOTO 83 

24850 CM0DE(NG2+III,I1,1)=VC31(I1) 

CMODE(NG2+III,I1,2)=O.D 00 
83 CONTINUE 

IF(NRGN.NE.l) GOTO 8612 
CRT (NG2+III, 1 )=XI 
NRT12=NRT/2 
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0119900 

0120000 

0120100 

0120200 

0120300 

0120400 

0120500 

0120S00 

0120700 

0120800 

0120900 

0121000 

0121100 

0121200 

1*121300 

0121400 

0121500 

0121600 

0121700 

0121800 

0121900 

0122000 

0122100 

0122200 

0122300 

0122400 

0122500 

0122600 

0122700 

0122800 

0122900 

0123000 

0123100 

0123200 

0123300 

0123400 

0123500 

0123600 

0123700 

0123800 

0123900 

0124000 

0124100 

0124200 

0124300 

0124400 

0124500 

0124600 

0124700 

0124800 

0124900 

0125000 

0125100 

0125200 


URITE(IN8,861 ) NG2,NRT12,III,CRT(NG2+III,1) 

URITE(IN8,87) <CM0DE(NG2+III,J,1), J=1,N) 

HRITE(IN8,87> (C*10DE(NG2+III, J,2) , J=1,N> 

8612 IF (NRGN.NE. 0 > GOTO 82 

861 F0RMAT(/3(1X,I4),1X,D15.8,1X,'REAL EIGEN VALUE AND MODE' ) 

WRITE (IN8, 861) NG2 , NRT , I r CRT(NG2+I,l> 

URITE(IN8,87) (CM0DE(NG2+I ,J,1),J=1,N) 

WRITE(IN8,87) CCM0DE(NG2+I ,0,2>,d=l,N) 

URITE(IN8,87) (CM0DECNG2+I+1, J, 1 > , J=1,N> 

WRITE (IN8, 87) (CM0DECNG2+I+1, J,2) , J=1,N) 

87 FORMAT(10(1X,D10.4)> 

82 CONTINUE 

IF(NRGN.EQ.l) NRT=NRT/2 
86 RETURN 
END 

FUNCTION HPX(N,I,P,A,B,YN,Y,MI,NRGN> 

IMPLICIT REAL (A-H,0-Z> 

DIMENSION P(6,4) 

CHO (X ) =DCOSH ( X ) * ( 1 -NRGN ) +DEXP ( A*X ) *NRGN 

SHO ( X ) =DSINH ( X ) * ( 1 -NRGN > +DEXP ( A*X ) *NRGN 

CHI (X)=SH0(X)*A 

SHI (X)=CH0(X>*A 

CH2(X)=CH0(X)*A*A 

SH2 ( X ) =SH0 ( X ) * A* A 

CH3 ( X ) =SH0 ( X ) *A*A* A 

SH3 ( X ) =CH0 ( X ) * A* A* A 

CH4(X)=CH0(X)*A*A*A*A 

SH4 ( X ) =SH0 ( X ) * A* A* A* A 

CH5 (X ) =SH0 (X )*A*A*A*A*A 

SH5 ( X ) =CH0 ( X > *A*A*A*A*A 

CH6 ( X ) =CH0 ( X > * A* A*A* A* A* A 

CH6 ( X > =SH0 ( X ) * A* A* A*A* A* A 

SHC 0 ( X ) =SH0 ( X ) *DCOS ( B*X ) 

SHS0(X)=SH0 (X)*DSIN(B*X) 

CHSO (X > =CH0 (X >*DSIN (B*X ) 

CHCO (X ) =CH0 ( X ) *DCOS ( B*X ) 

SHC1 ( X > =A*CHC0 ( X ) -B»SHS0 ( X ) 

SHS1 (X ) =A*CHS0 (X )+B*SHCO (X ) 

CHS1 (X ) =A*SHSO (X )+B*CHCO (X ) 

CHC1 (X > =A*SHC0 (X ) -B*CHS0 (X ) 

SHC2 (X ) = ( A*A-B*B )*SHC0 ( X ) -2*A*B»CHS0 (X ) 

SHS2 (X ) = ( A»A-B*B)*SHS0 (X >+2*A*B*CHC0 (X ) 
CHS2(X)=(A*A-B*B)*CHS0(X)+2*A*B*SHC0(X> 

CHC2 (X ) = ( A*A-B*B)*CHCO (X >-2*A»B*SHS0 (X ) 
SHC3(X)=(A*»3-3*A*B*B)*CHC0(X)-(3*A*A*B-B»*3)*SHS0<X) 

SHS3 (X ) = ( A**3-3*A*B*B)*CHS0 (X )+ (3*A*A*B-B**3 )*SHCO (X) 
CHS3(X)=(A**3-3*A*B*B)*SHS0(X)+(3*A*A*B-B**3)*CHC0(X) 

CHC3 (X ) = ( A**3-3*A*B*B ) *SHCO ( X ) - ( 3*A* A»B-B*»3 ) *CHSO ( X ) 
SHC4(X)=(A*<*4-6*A*A*B*B+B**4)*SHC0(X)-(4*A**3*B-4*A*B**3)*CHS0(X) 
SHS4 (X ) = ( A**4-6* A*A*B*B+B**4 )*SHSO (X )+ (4*A**3*B-4*A*B**3 >*CHCO ( X ) 
CHS4 (X ) = (A**4-6*A* A*B*B+B**4 )*CHSO (X )+ (4*A»*3*B-4*A«B**3 )*SHCO CX ) 
CHC4 (X ) = (A**4-6*A*A*B*B+B»*4)*CHC0 (X)- (4*A**3*B-4*A*B*»3)*SHS0 (X ) 
SHC5 (X ) = ( A**4-6* A* A*B*B+B**4 ) *SHC1 ( X ) - ( 4*A**3*B-4*A*B**3 ) *CHS1 ( X ) 
SHS5CX ) = ( A**4-6*A*A*B*B+B**4 )*SHS1 (X )+ (4*A**3*B-4*A*B**3 ) *CHC1 ( X ) 
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0125300 CHS5(X>=(A**4-6*A*A*B*B+B**4)*CHS1(X)+<4*A**3*B-4*A*B*»3)*SHC1(X) 

0125400 CHCS<X>=(A*<*4-6*A*A*B*B+B**4)*CHC1(X>-(4*A**3*B-4*A*B**3)*SHS1(X) 

0125500 SHC6 (X ) = ( A**4-6*A*A*B*B+B**4 >*SHC2 (X > - (4*A**3*B-4*A‘*B**3 >*CHS2 ( X ) 

0125600 SHS6 ( X ) = ( A**4-6»A* A*B*B+B**4 ) *SHS2 ( X ) + ( 4* A*»3*>B-4*A*B**3 ) *CHC2 ( X ) 

0125700 CHS6(X)=(A**4-6*A*A*B*B+B**4)*CHS2(X)+(4*A**3*B-4*A*B**3)*SHC2(X) 

0125800 CHC6 (X ) = ( a**4-6*A*A*B*B+B**4 )*CHC2(X ) - (4*A**3»B-4*A*B**3 >*SHS2 (X ) 

0125900 N1=N+1 

0126000 G0T0(1,2,3,4),N1 

0126100 1 IF(nr.EQ.l>HPX=(P(I,l)+P(I,2>*Y+P(I,3>*Y*Y+P(I,4)*Y*Y*Y> 

0126200 IF(MT.EQ.2)HPX=(P(I,1)*CH0(Y)+P(I,2)*SH0(Y) ) 

0126300 IF(HT.EQ.3>HPX=(P(I,1>*CHC0(Y)+P(I,2)*CHS0(Y)+P(I,3)*SHC0(Y)+P(I,4 

0126400 •)*SHS0(Y)> 

0126500 RETURN 

0126600 2 IF(MT.EQ.1)HPX=(P(I,2>*1 +P(I,3)*2*Y+P(I,4)*3*Y*Y)*YN 

0126700 IF (HT.EQ.2)HPX=(P(I,1)»CH1 (Y)+P(I,2)*SH1 (Y ) )*YN 

0126800 IF(MT.EQ.3)HPX=(P(I,1)*CHC1(Y)+P(I,2)*CHS1(Y)+P(I,3)*SHC1 (Y)+P(I,4 

0126900 *>*SHS1(Y>>*YN 

0127000 RETURN 

0127100.3 IF(MT.EQ.1>HPX=( P(I,3)*2 +P(I,4)*3*2*Y)*YN**2 

0127200 IF(MT.EQ.2)HPX=(P(I,1>*CH2<Y>+P(I,2)*SH2(Y)>*YN**2 

0127300 IF (MT.EQ.3)HPX=(P(I,1)*CHC2(Y)+P(I,2)*CHS2(Y)+P (I,3)*SHC2(Y)+P(I,4 

0127400 •)*SHS2(Y) )*YN**2 

0127500 RETURN 

0127600 4 IF(MT.EQ.1>HPX=( +P(I,4)*3*2*1)*YN**3 

0127700 IF (MT.EQ.2>HPX=(P(I,1)*CH3(Y)+P(I,2)*SH3(Y) )*YN**3 

0127800 IF(MT.EQ.3>HPX=(P(I,l>»CHC3(Y)+P(I,2)»CHS3(Y)+Pa,3)*SHC3(Y)+P(I,4 

0127900 •>*SHS3(Y)>*YN**3 

0128000 RETURN 

0128100 END 

0128200 SUBROUTINE FUN(NL,N,I1,N2,N4,NRGN,NV,MT,NMT,H,Y1,Z1 ,F,FPS, 

0128300 •CB11,CB12,CB22,SNC2,SHM2,SHQ2,KQ2,KM2,KP2,KN2,EPSL0N,CURVTR, THICK, 

0128400 *NHTY PE , SLN , CRT , CHODE , RCH , NCN1 , NCN2 , IPR ) 

0128500 IMPLICIT REAL (A-H.O-Z) 

0128600 REAL RCM 

0128700 REAL KM1,KM2,KN1 ,KN2,KQ1,KQ2,KP1,KP2,NU12,NU13,NU23 

0128800 DIMENSION CBll(l) ,CB12(1 ),CB22(1 ) ,SNC2(1> ,SHM2(1) ,SHQ2<1 ) ,KQ2(1 ), 

0128900 *KM2 ( 1 ) , KP2 ( 1 ) , KN2 ( 1 ) , THI CK ( 1 ) , NMTYPE ( 1 ) , 

0129000 *SLN(l),CRT(NCNl,l),CttODE(NCNl,NCN2,l),RCM(l),P(6,4) 

0129100 DATA Z01,Z02,Z03,Z04,Z05,Z06/0. ,0 . ,0. ,0. ,0. ,0./ 

0129200 C NV=1 N2 (XI, X2) NV=2 U2(X1,X2,Z) 

0129300 C NV=3 M2 <X1,X2) NV=4 PHI(X1,X2,0) 

0129400 C NV=5 Q2 (XI, X2) NV=6 H (X1,X2,Z) 

0129500 C NV=7 SYY(X1,X2,Z) NV=8 SYZ(X1,X2,Z) NV=9 SZZ(X1,X2,Z) 

0129600 C NV=10 HB(X1,X2> THICKNESS STRETCH PARAMETER 

0129700 C NV=11 SYZ(X1,X2,-1),X2 SHEAR STRESS DERIVATIVE AT PLY BOTTOM 

0129800 N1=NMTYPE(I1) 

0129900 CALL PZER0(P,24) 

0130000 CN=THICK(Il)/2. 

0130100 YN=CN/(H/2. ) 

0130200 Y=Yl/(H/2. ) 

0130300 Z=Z1/CN 

0130400 A=CRT(N2,1> 

0130500 B=CRT(N2,2) 

0130600 IF(4*N-5.LT.3.0R.4*N-5.GT.4*NL-2> GOTO 8101 
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0130700 

0130800 

0130900 8101 

0131000 

0131100 

0131200 

0131300 

0131400 8201 

0131500 

0131600 

0131700 8945 

0131800 9905 

0131900 

0132000 8946 

0132100 

0132200 

0132300 

0132400 8947 

0132500 9904 

0132600 1 

0132700 

0132800 

0132900 

0133000 

0133100 

0133200 

0133300 

0133400 3 

0133500 

0133600 

0133700 

0133800 

0133900 

0134000 

0134100 

0134200 5 

0134300 

0134400 

0134500 

0134600 

0134700 

0134800 

0134900 

0135000 

0135100 

0135200 

0135300 

0135400 

0135500 4 

0135600 

0135700 

0135800 

0135900 6 

0136000 


205=SLN(N4+4*N-5) 

Z06=SLN(N4+4*N-4) 

Z01=SLN(N4+4*N-3> 

Z02=SLN(N4+4*N-2) 

IF(4*N-l.LT.3.0R.4*N-l.GT.4*NL-2) GOTO 8201 
Z03=SLN(N4+4*N-1> 

Z04=SLN<N4+4*N-0> 

IF(4*N-5.LT.3.0R.4*N-5.GT.4*NL-2) GOTO 9905 

DO 8945 1=1,4 

P(5,I)=CH0DE(N2,4*N-5,I) 

P(6,I>=CM0DE(N2,4*N-4,I> 

DO 8946 1=1,4 
P(1,I)=CM0DE(N2,4*N-3,I> 

P(2,I)=CM0DE(N2,4*N-2,I) 

IF(4*N-l.LT.3.0R.4*N-l.GT.4*NL-2) GOTO 9904 

DO 8947 1=1,4 

P(3,I)=CM0DE(N2,4*N-1,I) 

P(4,I)=CW0DE(N2,4*N-0,I) 

G0T0(1,2,3,5,5,6,7,5,9,6,11),NV 

FPS=2*CN* ( KP2 (N1 ) * (Z06+Z04 ) +CB12 (N1 ) *EPSLON ) 

H11=HPX(1,1,P,A,B,YN,Y,HT,NRGN) 

H15=HPX(1,5,P,A,B,YN,Y,KT,NRGN> 

H13=HPX(1,3,P,A,B,YN,Y,HT,NRGN> 

H06=HPX(0,6,P,A,B,YN,Y,MT,NRGN> 

H04=HPX(0,4,P,A,B,YN,Y,HT,NRGN) 

F=2«CN* (H11*CB22(N1 )+ (H15-H13 )*KN2(N1 )+ (H06+H04 )*KP2 (N1 ) ) 

RETURN 

FPS= (KQ2 (N1 )/CN* (Z06-Z04 ) -CB12 (N1 )*CURVTR)*2 . /3*CN**3 
H22=HPX(2,2,P, A,B,YN,Y,HT,NRGN) 

H15=HPX(1,5,P,A,B,YN,Y,HT,NRGN) 

H13=HPX(1,3,P,A,B,YN,Y,HT,NRGN) 

H06=HPX(0,6,P,A,B,YN,Y,WT,NRGN) 

H04=HPX(0,4,P, A,B, YN,Y,tIT,NRGN) 

F=2 . /3*CN*CN* ( -H22*CB22 (N1 )+ (H15+H13 )*K«2 (N1 )+ (H06-H04 )*KQ2 ( N1 ) ) 
RETURN 

FPS=CN* (Z05+Z03 ) 

H32=HPX(3,2,P,A,B, YN,Y,HT,NRGN) 

H25=HPXC2,5,P,A,B,YN,Y,ttT,NRGN) 

H23=HPX (2, 3, P, A,B, YN, Y, HT, NRGN ) 

H16=HPX(l,6,P,A,B,YN,Y,ttT,NRGN) 

H14=HPX(1,4,P,A,B,YN,Y,HT,NRGN) 

H05=HPX(0,5,P,A,B,YN,Y,tIT,NRGN) 

H03=HPX(0,3,P,A,B,YN,Y,HT,NRGN) 

F=2 . /3*CN* ( -H32*CB22(N1 )+ (H25+H23 )*KM2 (N1 )+(H16-H14 )*KQ2 (N1 ) ) 

•+ ( H05+H03 ) *CN 
IF (NV.EQ.4) GOTO 4 
IF (NV.EQ.8) GOTO 8 
RETURN 

FPS=(SHH2(N1 >+SHQ2 (N1 ) )* (Z05+Z03) 

H12=HPX(l,2,P,A,B,YN,Y,!tT,NRGN) 

F=SHQ2 (N1 )/CN*F-H12+SHM2 (N1 )* (H05+H03) 

RETURN 

ZZZZ=.000000001D 00 
IF(NV.EQ.IO) Z=ZZZZ 
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0136100 

0136200 

0136300 

01364100 

0136500 

0136600 

0136700 

0136600 

0136900 

0137000 

0137100 

0137200 

0137300 

01374100 

0137500 

0137600 

0137700 

0137800 

0137900 

0138000 

0138100 

0138200 

0138300 

01384100 

0138500 

0138600 

0138700 

0138800 

0138900 

0139000 

0139100 

0139200 

0139300 

0139400 

0139500 

0139600 

0139700 

0139800 

0139900 

01410000 

01410100 

0140200 

0140300 

0140400 

0140500 

0140600 

0140700 

0140800 

0140900 

0141000 

0141100 

0141200 

0141300 

0141400 


61 FPS=WVD ( 1 , N1 , Z , NMT , RCM, IPR ) *CH*Z02+HVD (3 , N1 , Z*CN , NOT , RCM , IPR ) 
•+( (WVD(7 r Nl,Z,NMT,RCM,IPR)+UVD(8 r Nl,Z,NMT f RCM,IPR) )*Z06 

• + (UVD(7,N1,Z,NMT,RCM,IPR)-UVD(8,N1,Z, NttT, RCH, IPR) )*Z04)*CN 
H02=HPX(0,2,P,A,B,YN,Y,MT,NRGN> 

H22=HPX(2 t 2,P,A,B,YN,Y,HT,NRGN) 

H11=HPX(1,1,P,A,B,YN,Y,MT,NRGN) 

H15=HPX(1,5,P,A,B,YN,Y,HT,NRGN) 

H13=HPX(l,3,P r A,B,YN,Y,MT,NRGN) 

H06=HPX(0,6,P,A,B,YN,Y,MT,NRGN) 

H04=HPX(0,4,P,A,B, YN, Y,MT,NRGN) 

F=CN* ( H02*HVD ( 1 , N1 , Z , NMT , RCM, IPR ) +H22*MVD ( 2 , N1 , Z , NMT , RCH, IPR ) + 

• Hll*UVD(4,Nl,Z,NttT,RCM,IPR)+ 

• H15+(UVD(5,N1,Z,NMT,RCM, IPR)+UVD(6,N1,Z,NMT,RCM, IPR) )+ 

• H13*(-UVD(5,Nl,Z,NMT r RCM,IPR)+WVD<6,Nl,Z,NMT,RCM,IPR))+ 

• H06*(HVD(7,N1,Z,NMT,RCM,IPR)+WVD(8,N1,Z,NMT,RCM,IPR))+ 

• H04*(UVD(7,N1,Z,NMT,RCM,IPR)-UVD(8,N1,Z,NMT,RCM,IPR) ) ) 

IF (NV.EQ.6) RETURN 

IF (Z.NE.ZZZZ) GOTO 62 

FPS1=FPS 

F1=F 

z=-zzzz 

GOTO 61 

62 FPS=(FPSl-FPS)/2/CN/ZZZZ 
F=(Fl-F)/2/CN/ZZZZ 
RETURN 

1 - FPS=0 

H13=HPX(1,3,P,A,B,YN,Y,MT,NRGN) 

F=H13/CN 

RETURN 

2 FPS=Z01*CN*HVD(11,N1,Z,NMT,RCM,IPR)+Z05*CN*(HVD(13,N1,Z,NMT, 

•RCM,IPR)+WVD(15,N1,Z,NMT,RCM, IPR) )+ 

•Z03*CN*(-WVD(13,Nl,Z,NttT,RCM,IPR)+UVD(15,Nl,ZNMT,RCH,IPR)) 
H12=HPX(1,2,P,A,B,YN, Y,MT,NRGN) 

H32=HPX(3,2,P,A,B,YN,Y,HT,NRGN) 

H01=HPX(0,1,P,A,B,YN,Y,MT,NRGN) 

H21=HPX(2,1,P,A,B, YN,Y,HT,NRGN) 

H05=HPX(0,5»P, A,B,YN,Y,MT,NRGN) 

H25=HPX(2,5,P,A,B,YN,Y,MT,NRGN) 

H03=HPX(0,3,P,A,B,YN,Y,MT,NRGN) 

H23=HPX(2,3,P,A,B,YN,Y,MT,NRGN) 

H16=HPX(1,6,P, A,B,YN,Y,MT,NRGN) 

H14=HPX(1,4,P,A,B,YN,Y,MT,NRGN) 

F=(H12*WVD(9,N1,Z,NHT,RCM,IPR)+H32*HVD(10,N1,Z,NMT,RCM,IPR)+ 

• H01*UVD(11,N1,Z,NMT,RCM,IPR)+H21*WVD(12,N1,Z,NMT,RCM,IPR)+ 

• H05*( WVD(13,N1,Z,NMT,RCM,IPR)+HVD(15,N1,Z,NMT,RCM, IPR) )+ 

• H03*(-UVD<13,N1,Z,NMT,RCM,IPR)+HVD(15,N1,Z,NMT,RCM,IPR))+ 

• H25*( UVD(14,N1,Z,NMT,RCM,IPR)+UVD(16,N1,Z,NHT,RCM, IPR) )+ 

• H23*(-WVD(14,Ni,Z,NMT,RCM,IPR)+WVD(16,Nl,Z,NMT,RCH,IPR))+ 

• H16*( WVD(17,N1,Z,NMT,RCM,IPR)+UVD(18,N1,Z,NMT,RCM,IPR))+ 

• H14*( HVD<17,N1,Z,NMT,RCM,IPR)-UVD(18,N1,Z,NMT,RCM,IPR)))*CN 
RETURN 

8 FPS=Z05* (+Z/2+ . 5 )+Z03« ( -Z/2+ . 5 ) 

H05=HPX(0,5,P,A,B,YN,Y,HT,NRGN) 

H03=HPX(0,3,P,A,B,YN,Y,MT,NRGN) 
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0141500 F=H05* (Z/2+ . 5 >+H03* ( -Z/2+ .5 )+ ( 1-Z**2 )/2*F 

0141600 RETURN 

01411700 9 FPS=( .5+Z/2-(Z**3-Z)/4)*Z06+< .5-Z/2+(Z**3-Z)/4)*Z04 

01411800 H15=HPX(1 ,5,P,A,B, YN, Y,MT,NRGN) 

0141900 H13=HPX(1,3,P,A,B, YN,Y,MT,NRGN) 

0142000 H06=HPX (0,6,P,A,B,YN,Y,MT, NRGN ) 

0142100 H04 =HPX (0,4,P,A,B,YN,Y,MT, NRGN ) 

0142200 F=H06* ( . 5+Z/2- <Z**3-Z >/4 >+H04* ( . 5-Z/2+ (Z**3-Z >/4 )+ 

0142300 »H15* ( - (Z**2-l )/4- (Z**3-Z )/4 )+H13* ( <Z**2-1 )/4- (Z**3-Z>/4> 

0142400 RETURN 

0142500 7 FPS=Z06*(KP2(N1 >+Z*KQ2(Nl )+(Z**3- .6*Z)/4*SNC2(N1 ) )+Z04*(KP2(Nl )-Z* 

0142600 *KQ2<N1)-<Z**3-.6*Z)/4*SNC2<N1)>+CB12<N1)*<EPSL0N-CURVTR*CN*Z> 

0142700 H15=HPX(1,5,P,A,B,YN,Y,MT,NRGN) 

0142800 H13=HPX(1,3,P,A,B,YN,Y,MT,NRGN) 

0142900 H06=HPX(0,6,P,A,B,YN,Y,MT,NRGN) 

0143000 H04=HPX(0,4,P,A,B,YN,Y,MT,NRGN) 

0143100 H11=HPX(1,1,P,A,B,YN,Y,HT,NRGN) 

0143200 H22=HPX(2,2,P, A,B, YN,Y,MT,NRGN) 

0143300. F= Hll *CB22(N1)-H22 *CB22(N1)*Z 

0143400 *+H15 *( KN2(N1)+Z*KM2(N1)+SNC2<N1)»< Z**2-l./3+Z**3-.6*Z)/4) 

0143500 *+H13 •<-KN2(Nl)+Z*KM2<Nl)+SNC2(Nl)*(-Z**2+l./3+Z**3-.6*Z>/4> 

0143600 *+H06 •( KP2<N1)+Z*KQ2<N1)+.25*SNC2(N1)*<+(Z*»3-.6*Z>>> 

0143700 *+H04 •( KP2(N1)-Z*KQ2(N1)+.25«SNC2<N1)*<-(Z**3-.6*Z>)) 

0143800 RETURN 

0143900 END 

0144000 SUBROUTINE DEQMAT ( NSYH2 , NPLY , NEQ , NBC , NG1 , NG2 , MR, NMT , NSTREI,H, 

0144100 *CB22,KQ2,KM2,KP2,KN2, THICK, SLN,CRT,CMODE, 

0144200 *SK0 , SKI , SK2 , SK3 , SK4 , SX , VC45 , FSLN , ISGN , RCM , 

0144300 *NCN1, NCN2 , NSREG , NMTYPE , NG4 , N25L , IPR ) 

0144400 IMPLICIT REAL <A-H,0-Z) 

0144500 REAL RCM 

0144600 REAL KM1,KM2,KN1,KN2,KQ1,KQ2,KP1,KP2,NU12,NU13,NU23 

0144700 DIMENSION CB22(1),KQ2(1),KM2(1),KP2(1),KN2(1),THICK(1),SLN(1), 

0144800 *CRT (NCN1 , 1 ) , CMODE (NCN1 , NCN2, 1 ) , SX (NSTREI , 1 ) , 

0144900 *VC45 ( 1), FSLN (1),ISGN(1), NMTYPE <1>,SK0(N25L ,1), 

0145000 *SK1 (N25L ,1),SK2(N25L ,1),SK3(N25L ,1),SK4(N25L ,1),RCM(1) 

0145100 C U1,U2,H AND Z, ARE NONDIMENSIONALIZED H.R.T CN 

0145200 C XI, X2 ARE NONDIMENSIONALIZED U.R.T. LAMINATE SEMITHICKNESS, H/2 

0145300 C EPSLN=UNIFORM STRAIN IN Xl-DIRECTION, Ul,l 

0145400 C CURVTR=NONDIMENSIONALIZED CURVATURE IN Xl-DIRECTION, H,XX 

0145500 CALL PZERO (SK0,N25L*N25L) 

0145600 CALL PZERO ( SKI, N25L*N25L) 

0145700 CALL PZERO (SK2,N25L*N25L) 

0145800 CALL PZERO (SK3,N25L*N25L) 

0145900 CALL PZERO (SK4,N25L*N25L) 

0146000 CALL PZERO ( FSLN, N25L) 

0146100 C CREATING OF THE DIFFERENTIAL EQUATION OPERATOR MATRIX 

0146200 C SK4 *LM**4+ SK3 *LM**3+ SK2 *LM**2+ SKI *LM1+ SKO 

0146300 DO 101 N=1,NPLY 

0146400 L=NMTYPE ( N+NG1 ) 

0146500 YN=THICK(N+NGl)/2/(H/2) 

0146600 CN=THICK ( N+NG1 ) /2 

0146700 IF(N.LT.NPLY) L1=NMTYPE(N+NG1+1 ) 

0146800 IF^N.LT.NPLY) YNl=THlCK(N+NGl+l)/2/(H/2) 
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0136900 IF(N.LT.NPLY) CN1=THICK<N+NG1+1 >/2 

0137000 SK2 (3*N-3 F 3*N-3 ) =2 . *YN*YN*CB22 ( L ) 

0137100 IF (3*N-5 . GT . 2 . AND . 3*N-5 . LT . 3*NPLY-3 ) SK2(3*N-3,3*N-5)=2.*YN*YN*KN2- 

0137200 *(L) 

0137300 IF(3*N-1.GT.2.AND.3*N-1.LT.3*NPLY 3) SK2(3*N-3,3*N-1>=-2*YN*YN<*KN2- 

0137300 *(L) 

0137500 IF<3*N-5.GT.2.AND.3*N-5.LT.3*NPLY-3) SK0(3*N-3,3*N-5>=1 

0137600 IF(3*N-l.GT.2.AND.3*N-l.LT.3*NPLY-3> SK0(3*N-3,3*N-1)=-1. 

0137700 IF(3*N-3. GT.3. AND. 3*N-3.LT.3*NPLY-3> SK1<3*N-3,3*N-3)=2.*YN*KP2(L) 

0137800 IF (3*N .GT.3. AND. 3*N ,LT.3*NPLY-3) SK1<3*N-3,3»N)=2.*YN*KP2(L) 

0137900 SK3 (3*N-2, 3*N-2 ) =-2 . /3*CB22 (L >*YN*YN*YN*YN 

0138000 IF(3*N-5.GT.2.AND.3*N-5.LT.3»NPLY-3) SK3(3*N-2 r 3»N-5)=2./3*KH2(L)*- 

0138100 *YN*YN*YN 

0138200 IF(3*N-l.GT.2.AND.3*N-l.LT.3*NPLY-3) SK3(3*N-2,3*N-1 )=2./3*KH2(L)*- 

0138300 »YN*YN*YN 

0138300 IF (3*N-3 . GT . 3 . AND . 3*N-3 . LT . 3*NPLY-3 ) SK2(3*N-2,3*N-3>=2./3*KQ2(L)*- 

0138500 *YN*YN 

0138600 IF (3*N .GT.3. AND. 3*N .LT.3*NPLY-3) SK2(3*N-2,3*N )=-2./3*KQ2(L)- 

0138700 »*YN*YN 

0138800 IF(3*N-5.GT.2.AND.3*N-5.LT.3*NPLY-3) SKI <3*N-2,3*N-5)=YN 

0138900 IF(3*N-l.GT.2.AND.3*N-l.LT.3*NPLY-3) SK1(3«N-2,3*N-1)=YN 

0139000 IF(3*N-3. GT.3. AND. 3*N-3.LT.3*NPLY-3) SK0<3*N-2,3*N-3)=1 

0139100 IF (3*N . GT . 3 . AND . 3*N .LT.3*NPLY-3> SK0(3*N-2,3*N >=-l 

0139200 IF(N.EQ.NPLY) GOTO 101 

0139300 SK0(3*N-l,3*N-2)=WD(l,L,-l.,NOT,RCn,IPR)*CN 

0139300 SK2(3*N-l,3*N-2)=HVD<2,L,-l.,NHT,RCn,IPR)*YN*YN*CN 

0139500 SKO ( 3*N- 1 , 3*N+2 > =- WVD ( 1 , LI , 1 . 0 , NHT , RCM, IPR ) *CN1 

0139600 SK2 (3*N-1 ,3*N+2 > =-UVD(2 , LI , 1 . 0 , NWT, RCM, IPR)*YN1*YN1*CN1 

0139700 SK1(3*N-1,3*N-3)=HVD<3,L,-1.,NMT,RCH,IPR)*YN*CN 

0139800 SKI (3*N- 1 , 3*N+1 >=-UVD (3 , LI , 1 . 0 , NHT, RCtl, IPR >*YN1*CN1 

0139900 IF (3*N-5 . GT . 2 . AND . 3*N-5 . LT . 3*NPLY-3 ) SKI <3*N- 1 , 3*N-5 ) =WVD <5 , L , - 1 . , - 

0150000 •NMT , RCM , IPR ) • YN*CN+ HVD(6,L,-1.,NHT,RCH,IPR)*CN*YN 

0150100 IF(3*N-l.GT.2.AND.3*N-l.LT.3*NPLY-3) SK1(3*N-1,3*N-1)= 

0150200 •- HVD(5,L,-1.,NKT,RCH,IPR>*CN*YN- UVD(5,L1 f 1. »NKT r RCM F IPR)*CNl*YNl- 

0150300 *♦ HVD(6 r L,-l . r NMT,RCH»IPR)*CN*YN- WVD(6,L1,1. ,NHT,RCM,IPR)*CN1*YN1 

0150300 IF ( 3*N+3 . GT . 2 . AND . 3*N+3 . LT . 3*NPLY-3 ) SKl(3*N-l,3*N+3>= WVD(5,L1,1.- 

0150500 •0,NHT,RCH,IPR)*CN1*YN1- WVD(6,L1,1.0,NMT,RCM,IPR)*CN1*YN1 

0150600 IF (3*N-3 . GT . 3 . AND . 3*N-3 . LT . 3*NPLY-3 ) SK0(3*N-l,3»N-3>= WVD(7,L,-1.- 

0150700 *0,NMT,RCM,IPR>*CN+ HVD(8,L,-1. ,NHT,RCH,IPR)*CN 

0150800 IF (3*N .GT.3. AND. 3*N .LT.3*NPLY-3) SK0(3*N~1,3*N )= WVD(7,L,-1.- 

0150900 •0,NHT,RCW,IPR>*CN- HVD(8,L,-1. ,NMT,RCH,IPR)*CN- HVD(7,L1,1.0,NHT, - 

0151000 •RCil,IPR)*CNl- UVD(8,Ll,1.0,NHT,RCn,IPR)*CNl 

0151100 IF ( 3*N+3 .GT.3. AND . 3*N+3 . LT . 3*NPLY-3 ) SKO ( 3*N- 1 » 3*N+3 ) =- UVD(7,L1.1.- 

0151200 *0 , NHT , RCH, IPR)*CN1+ HVD(8,Ll,1.0,NnT,RCH,IPR)*CNl 

0151300 FSLN ( 3*N- 1 ) =- WVD(3,L,-CN,NHT,RCn, IPR)+UVD<3,Ll,CNl,NnT,RCH,IPR> 

0151300 SKl(3*N,3*N-2>= UVD(9,L,-1. , NWT, RCtt, IPR )*YN*CN 

0151500 SK3(3*N,3*N-2>= HVD(10,L,-1. ,N«T,RCH F IPR)*YN*YN*YN*CN 

0151600 SKI (3*N,3*W+2>=- WD<9,Ll,1.0,NHT,RCtt,IPR)*YNl*CNl 

0151700 SK3 ( 3*N , 3*N+2 ) =- HVD(10»Ll,1.0 F NnT,RCM F IPR)*YNl*YNl*YNl«CNl 

0151800 SK0(3*N,3*N-3>= HVD(ll F L r -l.»N«T»RCM,IPR)*CN 

0151900 SK2(3*N,3*N-3)= HVD(12,L r -l . ,NHT,RCM,IPR)*YN*YN*CN 

0152000 SKO (3*N, 3*N+1 )=- UVD(11»L1,1.0,NWT F RCH F IPR)*CN1 

0152100 SK2(3*N,3*N+1>=- HVD(12,L1,1.0,NHT F RCH F IPR)*YN1*YN1»CN1 

0152200 IF ( 3*N~5 . GT . 2 . AND . 3*N-5 . LT . 3*NPLY~3 ) SK0(3*N,3*N-5)= HVD(13,L,-1. 
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0152300 

0152400 

0152500 

0152600 

0152700 

0152800 

0152900 

0153000 

0153100 

0153200 

0153300 

0153400 

0153500 

0153600 

0153700 

0153800 

0153900 

0154000 

0154100 

0154200 

015430* 

0154400 

0154500 

0154600 

0154700 

0154800 

015490C 

0155000 

0155100 

0155200 

0155300 

0155400 

0155500 

0155600 

0155700 

0155800 

0155900 

0156000 

0156100 

0156200 

0156300 

0156400 

0156500 

0156600 

0156700 

0156800 

0156900 

0157000 

0157100 

0157200 

0157300 

0157400 

0157500 

0157600 


*NOT,RCH,IPR)*CN+ WVD(15,L,-1.,NOT,RCH,IPR)*CN 
IF(4*N-l.GT.2.AND.4*N-l.LT.4*NPLY-4) SK0(4*N,4*N-1)=- WVD(13,L,-1. 
*,NOT,RCH,IPR)*CN+ WVD(15,L,-1. ,NOT,RCH,IPR)*CN- WVD(13,L1,1.0,NOT, 
*RCN,IPR>*CN1-WVD(15,L1,1.0,NOT,RCH,IPR)*CN1 
IF ( 4*N+3 . GT . 2 . AND . 4*N+3 . LT . 4*NPLY -4 > SK0(4*N,4*N+3)= WVD(13,L1 , 1 . 0 
*,NOT,RCH,IPR)*CN1- WVD(15,L1,1.0,NOT,RCH,IPR)*CN1 
IF(4*N-5.GT.2.AND.4*N-5.LT.4*NPLY-4) SK2(4*N,4*N-5)=( WVD(14,L,-1. 
• , NUT , RCM , IPR )+WVD ( 16 , L , -1 . , NOT , RCI1 , IPR ) ) *YN* YN*CN 
IF(4*N~l.GT.2.AND.4*N-l.LT.4*NPLY-4) SK2(4*N,4*N-1 )= 

*(- WVD(14,L,-1. ,NOT,RCH,IPR)+ WVD(16,L,-1.,NOT,RCH,IPR))*YN*YN*CN- 
•(WVD(14,L1,1.0,NOT,RCH,IPR)+WVD(16,L1,1. , NOT, RCH, IPR) )*YN1*YN1*CN1 
IF (4*N+3 . GT . 2 . AND . 4*N+3 . LT . 4*NPLY-4 ) SK2(4*N,4*N+3)=( WVD(14,L1,1. 
*0 , NOT , RCM, IPR ) -UVD ( 16 , LI , 1 . 0 , NHT , RCH, IPR > ) *CN1*YN1*YN1 
IF ( 4*N-4 . GT . 3 . AND . 4*N-4 . LT . 4*NPLY-3 ) SKI (4*N,4*N-4)=( WVD(17,L,-1. 
* , NOT , RCtt , IPR ) +WVD ( 18 , L , - 1 . , NHT , RCH , IPR ) ) * YN*CN 
IF (4*N .GT.3.AND.4*N .LT.4*NPLY-3) SK1(4*N,4*N )= 

• ( WVD(17,L,-1. , NOT, RCH, IPR)- WVD(18,L,-1. ,NOT, RCtt, IPR) )*YN*CN- ( 
*WVD(17,Ll,1.0,NOT,RCn,IPR)+WVD(18,Ll,1.0,NOT,RCH,IPR))*YNl*CNl 
IF ( 4*N+4 . GT . 3 . AND . 4*N+4 . LT . 4*NPLY-3 > SKI (4*N , 4*N+4 ) =- ( UVD ( 17 , LI , 1 . 
•0, NOT, RCH, IPR)-WVD(18,L1 , 1 .0, NHT, RCH, IPR) )*YN1*CN1 
101 CONTINUE 
C NQQ=8 

C DO 6491 I=1,NQQ 

C6491 WRITE (6, 6390) I, (SK0(I, J), J=1,NQQ) 

C DO 6492 I=1,NQQ 

C6492 WRITE (6, 6390) I, (SK1(I,J),J=1,NQQ) 

C DO 6493 I=1,NQQ 

C6493 WRITE (6, 6390) I, (SK2(I, J) , J=1,NQQ) 

C DO 6494 I=1,NQQ 

C6494 WRITE (6, 6390) I, (SK3(I, J), J=1,NQQ) 

C DO 6495 1=1 ,NQQ 

C3495 WRITE (6, 6390) I, (SK4(I, J), J=1,NQQ) 

IF ( NSYH2 . GT . 1 . AND . NPLY . EQ . 1 ) GOTO 142 
IF(NSYH2.GT.1> GOTO 7104 
SLN(NG4+2)=0.D 00 
DO 5534 1=2, NPLY 
IF (NPLY.LT.2) GOTO 5534 

SLN (NG4+4*I-2 ) = (FSLN (4*1-5 )-SK0 (4*1-5, 4*I-6)*SLN(NG4+ 4*1-6))/ 

•SKO (4*1-5, 4*1-2) 

5534 CONTINUE 

DO 5536 1=1, NPLY 

5536 SLN(NG4+ 4*I-2)=SLN(NG4+ 4*I-2)/(THICK(NGl+I )/2. ) 

IF(NPLY.GT.l) GOTO 7113 
CHODE (NG2+5 , 1 , 1 ) =1 
CHODE ( NG2+2 , 2 , 1 ) =1 
CHODE (NG2+1 , 2 , 2 ) =1 
CHODE ( NG2+6 , 1 , 2 ) =1 
CH0DE(NG2+4,2,3)=1 
CHODE ( NG2+3 , 2 , 4 ) =1 
7113 IF(NPLY.LE.l) GOTO 142 
CHODE (NG2+5, 1,1) =1 
CHODE ( NG2+2 , 2 , 1 ) = 1 
CHODE (NG2+1 , 2, 2 ) =1 
CH0DE(NG2+6,1,2)=1 
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0157700 

0157800 

0157900 

0158000 

0158100 

0158200 

0158300 

0158400 

0158500 

0158600 

0158700 

0158800 

0158900 

0159000 

0159100 

0159200 

0159300 

0159400 

0159500' 

0159600 

0159700 

0159800 

0159900 

0160000 

0160100 

0160200 

0160300 

0160400 

0160500 

0160600 

0160700 

0160800 

0160900 

0161000 

0161100 

0161200 

0161300 

0161400 

0161500 

0161600 

0161700 

0161800 

0161900 

0162000 

0162100 

0162200 

0162300 

0162400 

0162500 

0162600 

0162700 

0162800 

0162900 

0163000 


C110DE(NG2+4,2,3)=1 
DO 7114 I=2,NPLY 
CllODE ( NG2+5 , 4*1 -3 , 1 ) =1 

CM0DE(NG2+2,4*I-2,1)=1 - 

CllODE(NG2+l,4*I-2,2)=-l./SKO (4*1-5, 4*I-2)*( 

* SKO (4*1-5, 4*1-6 )*C110DE(NG2+1, 4*1-6, 2) ) 
CllODE (NG2+1, 4*1-3, 1)=-1. /SKO (4*1-4, 4*1-3 )*( 

* SKO (4*1-4, 4*1-7 )*CM0DE(NG2+1, 4*1-7, 1 )+ 

* SKI (4*1-4, 4*I-6)*CM0DE(NG2+l,4*I-6,2)+ 

* SKI (4*1-4, 4*1-2 )*C110DE(NG2+1, 4*1-2, 2) ) 
CllODE (NG2+6, 4*1-3, 2)=-l. /SKO (4*1-4, 4*1-3 )*( 

* SKO (4*1-4, 4*1-7 )*CM0DE(NG2+6, 4*1-7, 2)) 
CllODE (NG2+6, 4*1-2, 1 )=-l. /SKO (4*1-5, 4*1-2 )*( 

* SKI (4*1-5, 4*1-7 )*CllODE(NG2+6, 4*1-7, 2 )+ 

* SKO(4*I-5,4*I-6)*CHODE(NG2+6,4*I-6,1)+ 

* SKI (4*1-5, 4*I-3)*CM0DE(NG2+6, 4*1-3, 2)) 
CM0DE(NG2+4, 4*1-2, 3)=-l./SK0(4*I-5,4*I-2)*( 

* SKO (4*1-5, 4*1-6) *CHODE ( NG2+4 , 4* I -6 , 3 ) ) 
CllODE (NG2+4, 4*1-3, 2)=-l. /SKO (4*1-4, 4*1-3 )*( 

* 2.*SKl(4*I-4,4*I-6)*CM0DE(NG2+4,4*I-6,3)+ 

* 2. *SK1 (4*1-4, 4*1-2 )*CM0DE(NG2+4, 4*1-2, 3 )+ 

* SKO (4*1-4, 4*I-7)*CM0DE(NG2+4, 4*1-7, 2) ) 
7114 CW0DE(NG2+4, 4*1-2, 1)=-1. /SKO (4*1-5, 4*1-2)* ( 

* SKO (4*1-5, 4*I-6)*CM0DE(NG2+4, 4*1-6, 1)+ 

* 2. *SK2 (4*1-5, 4*1-6 )*CH0DE(NG2+4, 4*1 -6,3)+ 

* 2.*SK2(4*I-5,4*I-2)*CW0DE(NG2+4, 4*1-2, 3)+ 

* SKI (4*1-5, 4*1-7 ) *CIi0DE ( NG2+4 , 4*1-7, 2)+ 

* SKI (4*1 -5 , 4*1 -3 )*CW0DE ( NG2+4 , 4*1 -3 , 2 ) ) 
CALL PZERO (SX,NSTREI) 

CALL PZERO (VC45,NSTREI) 

SX(1,1)=1 
SX(2,3)=1 
SX(3,4)=1 
VC45(3)=1 
DO 7116 I=1,NPLY 

SX (3+5*I-4, 5*1-3) =2. *SK2 (4*1-3, 4*1-3) 

IF (5*I-5.GT.4. AND. 5*1-5. LT.NSTREI+1) 

•SX (3+5*I-4, 5*1-5 )=SKO (4*1-3, 4*1-5) 

IF (5*1 .GT. 4. AND. 5*1 .LT.NSTREI+1) 

*SX(3+5*I-4,5*I )=SKO (4*1-3, 4*1-1) 

IF(I.EQ.NPLY) GOTO 7116 
SX (3+5*I-3, 5*1-1 )=SKO (4*1-1, 4*1-2) 

SX(3+5*I-3,5*I+4)=SK0 (4*1-1, 4*1+2) 

SX (3+5*I-2, 5*1-1 )=6.*SK2 (4*1-1, 4*1-2) 

SX(3+5*I-2, 5*1-2)= SKO (4*1-1, 4*1-2) 

SX(3+5*I-2,5*I-3)=2.*SK1(4*I-1, 4*1-3) 

SX (3+5*I-2, 5*1+4 )=6.*SK2 (4*1-1, 4*1+2) 

SX(3+5*I-2,5+I+3)= SK0(4*I-1, 4*1+2) 

SX(3+5*I-2,5*I+2)=2.»SK1 (4*1-1, 4*1+1) 

SX(3+5*I-1, 5*1-1 )=3.*SK1 (4*1-0, 4*1-2) 

SX(3+5*I-1, 5*1-3)= SK0(4*I-O, 4*1-3) 

SX(3+5*I-1,5*I+4)=3.*SK1 (4*1-0, 4*1+2) 

SX(3+5*I-1, 5*1+2)= SKO (4*1-0, 4*1+1) 

SX (3+5*I-0, 5*1-1 >=6.*SK3 (4*1-0, 4*1-2) 
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0163100 

0163200 

0163300 

0163400 

0163500 

0163600 

0163700 

0163800 

0163900 

0164000 

0164100 

0164200 

0164300 

0164400 

0164500 

0164600 

0164700 

0164800 

0164900 

0165000 

0165100 

0165200 

0165300 

0165400 

0165500 

0165600 

0165700 

0165800 

0165900 

0166000 

0166100 

0166200 

0166300 

0166400 

0166500 

0166600 

0166700 

0166800 

0166900 

0167000 

0167100 

0167200 

0167300 

0167400 

0167500 

0167600 

0167700 

0167800 

0167900 

0168000 

0168100 

0168200 

0168300 

0168400 


7116 


7117 

7104 


7119 


SX ( 3+5*1- 0, 5*1-2)= SKI (4*1-0, 4*1-2) 
SX(3+5*I-0,5*I-3)=2.*SK2(4*I-0, 4*1-3) 
SX(3+5*I-0, 5*1-4)= SKO (4*1-0, 4*1-3) 
SX ( 3+5*1- 0, 5*1+4) =6. *SK3 (4*1-0 4*1+2) 
SX(3+5*I-0, 5*1+3)= SKI (4*1-0 ,4*1+2) 
SX (3+5*1- 0, 5*1+2) =2. *SK2 (4*1-0, 4*1+1 ) 
SX(3+5*I-0, 5*1+1 )= SKO (4*1-0, 4*1+1 ) 
IF (5*I-5.GT.4. AND. 5*1-5. LT.NSTREI+1 ) 


•SX ( 3+5* I - 0 , 5*1-5)= SKO (4*1-0, 4*1-5) 

IF (5*1 .GT. 4. AND. 5*1 .LT.NSTREI+1) 

•SX(3+5*I-0,5*I )= SK0(4*I-0, 4*1-1) 

IF(5*I+5.GT. 4. AND. 5*1+5. LT.NSTREI+1) 

•SX(3+5*I-0, 5*1+5)= SKO (4*1-0, 4*1+3) 

CONTINUE 

CALL SKINV ( SX , NSTREI , SX , NSTREI , VC45 , NSTREI , 1 , SIMUB, NSTREI ) 
DO 7117 I=1,NPLY 

CM0DE(NG2+3,4*I-3,l)=VC45(5*I-4) 

CMODE (NG2+3,4*I-3,3) = VC45 (5*1-3) 

CMODE ( NG2+3 , 4* I -2 , 2 ) =VC45 ( 5*1 -2 ) 

CMODE ( NG2+3 , 4*1 -2 , 4 ) = VC45 (5* I - 1 ) 

IF(I.EQ.NPLY) GOTO 7117 

CMODE ( NG2+3 , 4*1 - 1 , 1 ) =VC45 (5* I ) 

CONTINUE 

IF (NSYM2.EQ.3) GOTO 142 
NEVN=(NPLY+1 )/2 


NEQ=4*NEVN-2 

IF (NPLY/2*2.EQ.NPLY) NEQ=4*NEVN 
N2EVN= (4*NPLY-2 ) /2+2 
N3EVN=(4*NPLY-2)/2+l 
N4EVN=(4*NPLY-2) 

IF ( NPLY/2*2 . NE . NPLY . AND . NSYM2 . EQ . 1 ) 
IF ( NPLY/2*2 . NE . NPLY . AND . NSYM2 . EQ . 1 ) 
IF ( NPLY/2*2 . NE . NPLY . AND . NS YM2 . EQ . 2 ) 
IF ( NPLY/2*2 . NE . NPLY . AND . NSYM2 . EQ . 2 ) 
IF ( NPLY/2*2 . EQ . NPLY . AND . NSYM2 . EQ . 1 ) 
IF ( NPLY/2*2 . EQ . NPLY . AND . NSYM2 . EQ . 1 ) 
IF ( NPLY/2*2 . EQ . NPLY . AND . NSYM2 . EQ . 2 ) 
IF ( NPLY/2*2 . EQ . NPLY . AND . NS YM2 . EQ . 2 ) 
IF ( NPLY/2*2 . NE . NPLY . AND . NSYM2 . EQ . 1 ) 
IF ( NPLY/2*2 . EQ . NPLY . AND . NSYH2 . EQ . 1 ) 
IF ( NPLY/2*2 . NE . NPLY . AND . NSYM2 . EQ . 2 ) 
IF ( NPLY/2*2 . EQ . NPLY . AND . NSYM2 . EQ . 2 ) 
DO 7119 I=1,NEVN 
0=4* (NPLY+l-I ) 

ISGN(J-3)= (4*I-3)*(-l.)**(NSYM2) 
ISGN(J-2)= (4*I-2)*(-l.)**(NSYM2+l) 
ISGN(J-5)= <4*I-1)*(-1.)**(NSYM2+1> 
ISGN(J-4)= (4*1-0 )*(-l.)**(NSYM2> 
DO 7121 I=1,N3EVN 
IF(I.EQ.I2) GOTO 7121 
DO 71212 0=N2EVN , N4EVN 
K=ISGN(J) 

JSGN=K/IABS(K) 

K=IABS(K) 


I1=N3EVN-1 

I2=N3EVN-1 

I1=N3EVN 

I2=N3EVN 

I1=N3EVN 

I2=N3EVN-1 

I1=N3EVN-1 

I2=N3EVN 

NBC = 8*(NEVN-l)+l+4 
NBC = 8* (NEVN-1 )+l+4+2+2 
NBC = 8* (NEVN-1 )+l+2 
NBC = 8* (NEVN-1 )+l+4+2 
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0168500 

0168600 

0168700 

0168800 

0168900 

0169000 

0169100 

0169200 

0169300 

0169400 

0169500 

0169600 

0169700 

0169800 

0169900 

0170000 

0170100 

0170200 

0170300- 

0170400 

0170500 

0170600 

0170700 

0170800 

0170900 

0171000 

0171100 

0171200 

0171300 

0171400 

0171500 

0171600 

0171700 

0171800 

0171900 

0172000 

0172100 

0172200 

0172300 

0172400 

0172500 

0172600 

0172700 

0172800 

0172900 

0173000 

0173100 

0173200 

0173300 

0173400 

0173500 

0173600 

0173700 

0173800 


SK0<I,K)=SK0<I,K)+SK0(I,J)*JSGN 
SK1(I,K)=SK1(I,K)+SK1<I,J)»USGN 
SK2(I,K)=SK2<I,K)+SK2(I,J)*JSGN 
SK3(I,K)=SK3(I,K)+SK3(I, J)*JSGN 
71212 SK4 ( I , K ) =SK4 ( I , K ) +SK4 ( I , J ) * JSGN 
7121 CONTINUE 

DO 7129 J=1,N3EVN 
SK0(I2,O)=0.D 00 
SK1(I2,J)=0.D 00 
SK2(I2,O)=0.D 00 
SK3(I2,0)=O.D 00 
7129 SK4(I2,J)=0.D 00 

CALL PZERO(SKO(1,I1),N3EVN) 

CALL PZERO ( SKI (1,11), N3EVN ) 

CALL PZERO (SK2(1, II ) ,N3EVN) 

CALL PZER0(SK3 (1,11), N3EVN ) 

CALL PZERO <SK4(1, II >,N3EVN) 

SK1(I2,I1)=1. 

6390 F0RMAT(1X,I2,14F9.4> 

C NQQ=8 

C DO 6391 1=1, NQQ 

C6391 WRITE (6, 6390) I, <SK0(I, J), J=1,NQQ) 

C DO 6392 1=1, NQQ 

C6392 WRITE (6, 6390) I, (SKI (I, J) , J=1,NQQ) 

C DO 6393 1=1, NQQ 

C6393 WRITE (6, 6390) I, (SK2(I,J),0=1,NQQ) 

C DO 6394 1=1, NQQ 

C6394 WRITE(6,6390 ) I, (SK3(I, J), J=1,NQQ) 

C DO 6395 1=1, NQQ 

C6395 WRITE(6,6390) I, (SK4(I, J) , J=1,NQQ> 

RETURN 

142 NEQ=4*NPLY-2 
NBC=8*NPLY-2 
RETURN 
END 

SUBROUTINE PZERO (A, N) 

IMPLICIT REAL (A-H.O-Z) 

DIMENSION A(l) 

DO 1 1=1, N 
1 A(I)=0.D 00 

RETURN 
END 

SUBROUTINE BID0T(A,I1,B,I2,C,I3) 

IMPLICIT REAL <A-H,0-Z) 

DIMENSION A(I1,1),B(I2,1),C(I1,1) 

DO 1 1=1,11 
DO 1 K=1,I3 
C(I,K)=0.D 00 
DO 1 0=1,12 

1 C(I,K)=C(I,K)+A(I,J)*B(J,K) 

RETURN 

END 

SUBROUTINE SKINV ( SK, NSK, SKI , NSKI , X , NX , INDIC , SIMUL , N ) 
IMPLICIT REAL (A-H,0-Z> 
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0173900 

C. . . 

. .1NDIC =0 GIVES DETE. (SIMUL), INVERSE (SKI) AND SOLUTION(X) 

0174000 

C. . . 

. .INDIC -VE GIVES DETERMINANT AND INVERSE 

0174100 

C • • • 

..INDIC +VE GIVES DETERMINANT AND SOLUTION 

0174200 


DIMENSION IRON (251 ) , JC0LC251 ) , J0RD(251 ) ,Y (251 ) , A(251 ,251 ) ,SK(NSK, 1 

0174300 


* ) , SKI ( NSKI ,1),X(NX> 

0174400 


EPS=l.D-70 

0174500 


NRC=251 

0174600 


MAX=N 

0174700 


IF ( INDIC. GE.O) MAX=N+1 

0174800 


DO 101 1=1, N 

0174900 


DO 101 J=1,N 

0175000 

101 

Ad, J)=SK(I, J) 

0175100 


IF(INDIC.LT.O) GO TO 103 

0175200 


DO 102 1=1, N 

0175300 

102 

A(I,MAX)=X(I> 

0175400 

103 

CONTINUE 

0175500 


IF (N.LT.NRC) GOTO 5 

0175600 


HRITE(6,200) 

0175700 

200 

FORMAT ( /IX,' MATRIX IS TOO BIG IN "SKINV" '/> 

0175800 


RETURN 

0175900 

5 

CONTINUE 

0176000 


DETER=1. 

0176100 


DO 18 K=1,N 

017S200 


KM1=K-1 

0176300 


PIVOT=0. 

0 ,-T S400 


DO 11 1=1, N 

0176500 


DO 11 J=1,N 

0176600 


IF(K.EQ.l) GOTO 9 

0176700 


DO 8 ISCAN=1,KM1 

0176800 


DO 8 JSCAN=1,KM1 

0176900 


IFd.EQ.IROUdSCAN) ) GOTO 11 

0177000 


IF(J.EQ.JCOL(JSCAN>> GOTO 11 

0177100 

8 

CONTINUE 

0177200 

9 

IF (DABS(A(I, J) ) .LE. DABS (PIVOT) ) GOTO 11 

0177300 


PIVOT=A(I,J) 

0177400 


IROW(K)=I 

0177500 


JCOL(K)=J 

0177600 

11 

CONTINUE 

0177700 


IF (DABS (PIVOT J.GT.EPS) GOTO 13 

0177800 


SIMUL=0 . 

0177900 


URITE(6,1) 

0178000 

1 

FORMAT (IX, 'STOPPED IN SUBROUTINE SKINV DETERMINANT = 0') 

0178100 


RETURN 

0178200 

13 

IROWK=IROU(K) 

0178300 


JCOLK=JCOL(K> 

0178400 


DETER=DETER*PIVOT 

0178500 


DO 14 J=1,MAX 

0178600 

14 

A (IROUK, J ) =A ( IROUK , J ) /PI VOT 

0178700 


A ( IROWK, OCOLK ) =1 . /PIVOT 

0178800 


DO 18 1=1, N 

0178900 


AIJCK=A(I,JCOLK) 

0179000 


IF(I.EQ. IROUK) GOTO 18 

0179100 


Ad, JCOLK) =- AI JCK/PIVOT 

0179200 


DO 17 J=1,HAX 
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0179300 

17 

IF(J.NE.JCOLK) Ad, J)=Ad, J)-AIJCK*AdROUK, J) 

0179400 

18 

CONTINUE 

0179500 


DO 20 1=1, N 

0179600 


IROUI=IROUd) 

0179700 


JCOLI=JCOLd> 

0179800 


JORD ( IROHI ) =JCOLI 

0179900 

20 

IF(INDIC.GE.O) X(JCOLI)=A(IROUI,HAX) 

0180000 


INTCH=0 

0180100 


NM1=N-1 

0180200 


DO 22 1=1, Ntll 

0180300 


IP1=I+1 

0180400 


DO 22 J=IP1,N 

0180500 


IF ( J0RD( J) .GE. JORD(I > ) GOTO 22 

0180600 


JTEMP=JORD(J> 

0180700 


JORD(J)=JORD(I) 

0180800 


JORD ( I ) =JTEMP 

0180900 


INTCH=INTCH+1 

0181000 

22 

CONTINUE 

0181100. 


IF(INTCH/2»2.NE.INTCH) DETER=-DETER 

0181200 


IF(INDIC.LE.O) GOTO 26 

0181300 


SIHUL=DETER 

0181400 


RETURN 

0181500 

26 

DO 28 J=1,N 

0181600 


DO 27 1=1, N 

0181700 


IROUI=IROUd> 

0181800 


JCOLI=JCOL(I) 

0181900 

27 

Y ( JCOLI ) =A (IROUI , J ) 

0182000 


DO 28 1=1, N 

0182100 

28 

Ad, J)=Y(I) 

0182200 


DO 30 1=1, N 

0182300 


DO 29 J=1,N 

0182400 


IROUJ=IROH(J> 

0182500 


JCOLJ=JCOL( J ) 

0182600 

29 

Y (IROUJ ) =A (I , JCOL J ) 

0182700 


DO 30 J=1,N 

0182800 

30 

Ad, J)=Y(J> 

0182900 


SIIIUL=DETER 

0183000 


DO 209 1=1, N 

0183100 


DO 209 J=1,N 

0183200 

209 

SKId, J)=A(I,J) 

0183300 


RETURN 

0183400 


END 
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APPENDIX II: PLY EQUATIONS 


The derivation of these equations Is reported In reference 8. These equa- 
tions encompass the following main categories: overall equilibrium, constitu- 
tive relations, stress distribution, and displacement distribution. 

Overall Equilibrium: 



,.k k 0 

N 2 , 2 * n 2 = 0 


(I-D 

H 2,2 + C k m 2 ~ Q 2 = 

0 

d-2) 

n k k ~ 

«2.2 * " ' 0 


d-3) 

where 



n 2 = °2z (x 2’ c k ) ' 

°2z< x 2’ - c k> 

(1-4) 

m 2 = °2z* x 2' C k ) 

a 2z ( V - c k> 

d-5) 

q = ° zz ( x 2 ' c k ) 

°zz^ x 2’ _c k ) 

(1-6) 


where N? f M 2 and Q 2 are, respectively, the force, moment, and shear 
resultants associated with the X2-coord1nate direction, and n2, m 2 , and 
q are related to the Interfaclal stresses. The semithickness of the k-th ply 
Is C|(. Superscript k which Identifies the generic ply will be dropped In 
the subsequent equations for convenience. 

Constitutive Relations: 

Nfl/h = C Ri e t C B2 U 2> 2 + K n |} c n 2j 2 + K p(3 P ; B = 1.2 (1-7) 

M li /I * -^I31 w ,ll - Ct32 u ,22 + K m0 m 2,2 f K qG q/c 0~ 8 ) 

<t>2 + w ,2 = (Q 2 - c m 2 /3) S 44 c2/(2I) (1-9) 

where 

K m 0 = (3 Cfl 2 S 3 a C a 2/C22 - 2 C 32 s 44 + 2Cg)/20 ; a, 6 = 1,2 
Kq|3 = O C 62 S 3a C a 2 /C 22 - 12 C 32 S 44 + 12C B )/20 

K n(i = (^B2 s 3a ^ q . 2^22 + C&2 S 44 + 2Cg)/12 
K p 0 Cp/2 

p = cf zz (x 2 , c) + o zz (x 2 , -c) , h = 2c , I = 2 c 3 /3 (1-10) 

where e Is the applied strain In the x-] - direction and W n Is the 
curvature about X2-ax1s. [C] Is the plane stress stiffness matrix and [S] 
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Is the general flexibility matrix. The constants C] , C?, and C B Introduce 
the Influence of the transverse normal stress In the computation on the refined 
surface parallel stresses. For a given lamina with arbitrary orientation of 
the fibers and 02 are surface parallel normal stresses, 03 Is 

transverse normal stress and 05 the In-plane shear stress. The corresponding 
engineering strains are denoted by e] , e 2, C3, and T5. Using the 
generalized Hooke's law, the equations for c \ , e 2 and T B can be 
Inverted to obtain 

0^ = *- 0^03 ; 1,j = 1, 2, and 6 ( I - 1 1 ) 

where_ e 5 Is an alternative notation for T B . The expressions for C] , C2, 
and C5 can be obtained In terms of flexibilities and later In terms of lamina 
elastic moduli with relative ease. 


The axial displacement U2 and the transverse displacement W are asso- 
ciated with reference surface z = 0. Here $2 Is the rotation of a normal to 
the reference surface. 

Stress Distribution: 


0 B = N B /h + M 0 z/I + K B n 2>2 (z 2 - c 2 /3)/(2h) 

r K b (c m 2> 2 + Q)(z 3 - 3 c 2 z/5)/(6I) 
0 2z = n 2 z/I * c m 2 ( 3 z 2 - c 2 )/(6I) Q 2 (z 2 - c 2 )/(2I) 

°zz = P/2 ~ n 2 , ?A z2 - c 2 )/(2h) + q z/h - (c (1^2 * q)(z 3 - c 2 z)/(6I) 
where 


K b = C B 2S 3a C a2 /C22 + C B 2 S 44 " c |3 i = 1 and 2 

01 = surface parallel stress In x-| direction, 

02 = surface parallel stress In X2 direction, 

©2 Z = shear stress In X2-Z plane, and 

a 7l = transverse normal stress In z direction. 
Displacement Distribution: 

w = W t S 3a (N a z/h t M a z 2 / ( 21 ) ) + S 3a K a n 2>2 (z 3 - c 2 z)/(6h) 
+ S 3a K a (c m 2> 2 * Q)(z 4 /4 - 3 c 2 z 2 /l 0)/( 61) 

+ S33(p z/2 <- q z 2 /(2h)) - S33n2 >2 (z 3 " 2 c 2 z)/(6h) 
r S 33 (c m 2f2 «- q)( z 4 /4 - c 2 z 2 /2)/(6I) 


( 1 - 12 ) 

(1-13) 

(1-14) 


(1-15) 


( 1 - 16 ) 
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U = u 2 - z W f2 - S 3 a (N a z2/(2h) + M a z 3 /( 6 I )) >2 

- S 3 aM"2,22 U4 / 4 ' c2z2/2 )/( 6h ) 

+ (c m 222 f q , 2 ) ( z5 /20 - c 2 z 3 /10)/(6I)| + S33 ^ p 2 z2/4 + Q,2 z3/ ( 6h ) 

- n 2 > 22 (z 4 /4 ~ 3 c 2 z 2 /2)/(6h) - (c m 222 + q > 2 )(zV20 - c 2 z 3 / 6 )/ 6 I)| 

* S 44 |z 2 n 2 / ( 2h) + c m 2 ( z 3 - c 2 z)/(61) - Q 2 (z 3 - 3 c 2 z)/( 6 I)| (1-17) 

The Greek subscripts In the foregoing expressions assume values of 1 
and 2. The repeated Greek Indlcles Indicate a summation over the appropriate 
range. Finally, u 2 and w are the displacement components In the x 2 
and z coordinate directions. 
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APPENDIX III: INPUT AND OUTPUT INSTRUCTIONS FOR THE CODE 


This code contalna a MAIN program and T 6 subroutines. Additionally, It 
also needs a subroutine to fine the determinant and Inverse of a complex 
matrix. The 1MSL mathematical subroutine library Is used for this purpose. 

The program Is restricted to cylindrical bending problems In X 2 -Z plane. 

The laminate can be analyzed In stretching In x-|-d1 rectlon or In bending In 
x-|-z plane. 

The MAIN program controls Input. The subroutine MPR 6 controls construc- 
tion of eigen solutions for each region, construction and solution of boundary/ 
continuity conditions matrix, and finally computation of response variables 
through the thickness and along the axis of the laminate. The core memory can 
be controlled by the dimension of RCM vector In MAIN program. Also, set the 
value of the variable MAXRCM to the same value as the dimension of RCM vector. 

The Input Is divided Into three parts. The first part given by Data cards 
1 through 5 describes the sublaminates, ply layup, and material properties. 

The second part generates the eigenvalues and eigen modes. Finally, boundary/ 
continuity conditions are Input through Data card 6 . The distribution of the 
response variables Is controlled by Data card 7. The first part of the Input 
Is explained below: 

Data 1 : NSREG, NMT 

Format( 12) .. .These variables, respectively, Indicate the number of sub- 
laminates which the laminate Is divided Into and the number of materials used. 

A different ply angle indicates a different material. 

Data 2: EPSLON, CURVTR 

Format(8E10.4) . . .EPSLON Is the uniform strain In x-| -dl rectlon and 
CURV1R Is the quotient of curvature In x-| - z plane. 

Data 3: NL, NGN, NGT 

Format(2014) . . .NL Is the number of layers Independently analyzed In a sub- 
laminate. A "0" for NGN constructs full solutions. A "1" gives decay type 
solutions for the current sublaminate. If NGT is either "0" or the current 
sublaminate Identification number, then the full eigenvalue problem Is solved 
for the current sublaminate. If NGT Is Initiated to a number less than the 
current sublaminate Identification number, then Its eigen solutions are not 
computed. Instead these are assumed to be the same as those of the NGT sub- 
laminate. This means that the present sublaminate Is Identically the same as 
the NGT sublaminate. The data for NSREG sublaminates Is given here. 

Data 4: (THICK(K), K=1,NL) 

Format(8E10. 4 ) .. .Thicknesses for the NL layers of the sublaminate. The 
data for all the NSREG sublaminates are entered here. 

Data 5: Eq , E?, E 3 , v\ 2 > ” 13 , 1 * 23 , G 23 , G 13 , Gq 2 » 6 

Format(10f8. 4) .. .Material properties along ply principal directions. 6 Is 
the angle of orientation. The data for all the NM1 materials Is entered here. 
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The second part of the data, to compute eigen modes. Is to be given Inter- 
actively. This procedure can be automated If a reasonable approximation of the 
eigen values Is known. It begins with the following series of prompts: 

Prompt 1: "GIVE 0 IF YOU 00 N01 WAN! P0LYN0M1NAL APPROX. 

NR = xx, nsym - xxx" 


A default ("0") uses complete eigen matrix to evaluate real roots. Any 
other number uses a polynomlnal representation. Control goes to Prompt 2. NR 
Indicates the current sublamlnate for which real roots are sought. There Is a 
do loop on NR covering all sublaminates. Sublamlnates are skipped for which 
NGN Is not Initiated to "0" or their proper Identification number. There Is a 
do loop on NSYM. Pure stretching modes are constructed If NGN Is 1. Pure 
bending modes are computed when It Is 2. Both are computed when NSYM = 3. 

These do loops end when the reply to Prompt 5 Is a default. 

Prompt 2: "GIVE 0 IF YOU WAN! REAL ROOTS" 

A default computes real roots by going to Prompt 9. Any other number 
transfers to Prompt 3. 

Prompt 3: "GIVE 0 IF YOU DO N01 WAN! POLYNOMIAL APPROX. 

NR = xx, nsym = xxx" 

A default ("0") used complete eigen matrix to evaluate complex roots. Any 
other number uses a polynomial representation. Control goes to Prompt 4. 

Prompt 4: "GIVE 0 IF YOU WANT COMPLEX ROOTS" 

A default computes complex roots by going to Prompt 13. Any other number 
transfers to Prompt 5. 

Prompt 5: "GIVE 0 IF POLYNOMIAL IS NOT DESIRED" 

A default ends do loops on NR and NSYM. Any other number Initiates the 
process for construction of a polynomial suppressing the real and complex roots 
constructed thus far by transferring control to Prompt 6. Here, It also shows 
the roots computed thus far. 

Prompt 6: "GIVE NRT FOR POLYNOMIAL. . .DEFAULT IS AS IS A NEGATIVE VALUE 

INITIATES NRT TO 0" 

NR1 Indicates the number of real roots computed thus far for the particu- 
lar cases of NR and NSYM. A negative value Initiates NRT to 0. Any other 
number less than the current NR1 Initiates NRT to that value. This helps con- 
trolling the number of real roots In the polynomial evaluation. Control next 
goes to Prompt 7. 

Prompt 7: "GIVE NCI FOR POLYNOMIAL. . .DEFAULT IS AS IS A NEGATIVE VALUE 

INIUA1ES NCI TO 0" 

NCI Indicates the number of complex roots computed thus far for the par 
tlcular cases of NR and NSYM. A negative value Initiates NCI to 0. Any other 
number less than the current NCI Initiates NCI to that value. This helps 
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controlling the number of complex roots In the polynomial evaluation. Control 
goes to Prompt 8. 

Prompt 8: "GIVE 0 IF POLYNOMIAL IS NOT DESIRED" 

A default transfers the control to Prompt 1 without changing either NR or 
NSYM. Any other number computes the polynomial. Control again shifts to 
Prompt 1 . 

Prompt 9: "10F 10.6. . .XI ... INITIAL ROOT VALUE" 

This commences the process for finding real roots. This Is based on 
Newton- Raphson technique. XI Is the Initial root value. A default transfers 
control to Prompt 11. If there Is convergence, the root Is displayed and the 
control goes to Prompt 10. If there Is no convergence after 100 Iterations the 
control goes to Prompt 9 displaying the fact. 

Prompt 10: "1 2 . . . INPUT . . . J=1 TO STORE ROOT" 

After the converged root Is displayed, It asks whether the root Is to be 
saved (Input "1") or not (default) for the purpose of a later storage. If 
allowed, It stores two values for the root, +xxx and -xxx, and appends NRT by 
2. If NGN Is not Initiated to "0" only one root Is stored and NRT Is Increased 
by 1. The choice X2~or1g1n In relation to the edge where decay type solu- 
tions are required determines whether to store positive or negative roots. 
Control goes to Prompt 9. NRT starts at 0 for each NR and NSYM. At the behest 
of Prompt 6 this parameter can be changed to drop unwanted roots. 

Prompt 11: "12. . .INPUT. . ,J=0 TO FIND REAL EIGEN MODES" 

A default computes the eigen modes. Next Is Prompt 12. 

Prompt 12: "GIVE 0 IF REAL ROOTS ARE TO BE STORED" 

A default stores real roots and real modes and sets NRT back to zero. 

Once the roots are stored, then the polynomial at Prompt 5 cannot be con 
structed by suppressing these roots. Control goes to Prompt 3. 

Prompts 13 through 16: 

The "In" control comes from Prompt 4 and the "out" control goes to 
Prompt 5. These prompts do the same job as described for Prompts 9 through 12, 
but for complex roots. NCT keeps track of the number of complex roots to be 
stored. Four roots. . .positive, positive complex conjugate, negative, and nega- 
tive complex conjugate versions will be saved and NCT will be appended by 4 If 
NGN Is "0." Otherwise, only positive or negative roots (depends on the discre- 
tion of the user) will be saved and NCT will be padded up by two. 

The logic outlined by these prompts Is elementary and Is useful for find- 
ing the eigen solutions. User ends this process after determining that all 
roots are found and stored. Once a root Is stored at the Instance of either 
Prompt 10 or 14, It will be suppressed so that the next Iteration will not con- 
verge on to this root. The next set of data concerns the boundary/continuity 
conditions and printing of response variables. 
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Data 6: NBCREG, NBCPLY, NBCV, YNBC , ZNBC, PNBC, NBCREGL 

Format( 314 , 3E1 2 . 6 , 14) . . . Repeat this data to cover all boundary and conti- 
nuity conditions. NBCREG Is the sublamlnate, NBCPLY Is the ply, NBCV Is the 
variable, YNBC Is the local X2-coord1nate ZNBC Is the local z coordinate. 
The origin for the local z axis Is always In the local ply middle plane. 
Origin for x 2 - coordl nate can be placed anywhere but care should be taken 
when using pure decay type solutions. For each sublamlnate, the plys are num 
bered 1, 2, 3,. ..locally from top to bottom. 


NBCV 


1 

Indicates 

N 2 (x 2 ) 

NBCV 


2 

Indicates 

u 2 (x 2 , z) 

NBCV 


3 

Indicates 

M 2 (x 2 ) 

NBCV 

- 

4 

Indicates 

4>2(x 2 ) 

NBCV 


5 

Indicates 

Q2( x 2) 

NBCV 

= 

6 

Indicates 

w(x 2 ,z) 

NBCV 

= 

7 

Indicates 

°22( x 2 * z ) 

NBCV 


8 

Indicates 

c 2z( x 2 > z ) 

NBCV 

= 

9 

Indicates 

°zz(*2» z ) 

NBCV 


10 

Indicates 

transverse stretch parameter 

NBCV 

= 

11 

Indicates 

°2z,2( x 2» z ) 


The tenth variable, transverse thickness stretch parameter, denotes the 
coefficient of z In the thlcknesswl se distribution of the transverse dis- 
placement. N 2 , M 2 , and Q 2 are defined layerwlse when full solutions are 
used. At edges where decay type solutions are present, this definition Is 
again retained. At far edges, however, these are defined over entire sublaml- 
nate thickness preserving local force equilibrium. For conditions Involving 
several variables (NV's), this data continues on to the next card and Is termi- 
nated by specifying NBCREGL = 0 on the last card. 

For example, 


2 

1 

1 . 1 OOOOOd 

01 

. OOOOOOd 

00 . 1 20000d 

03 

8 

1 

3 

3- . 1 OOOOOd 

01 

. OOOOOOd 

00-.1 OOOOOd 

01 

8 

2 

2 

2 . 200000d 

01- 

. 200000d 

00 .900000d 

01 

8 

4 

7 

5-. 1 OOOOOd 

01- 

. 200000d 

00 . 200000d 

01 

0 


Indicates the following: 

N 2 - M 2 * 9U 2 + 2Q 2 + 120 = 0 
where, 

N 2 Is defined for the second sublamlnate and first ply at X 2 = 1 , 

M 2 Is defined for the first sublamlnate and third ply at X 2 = -1 , 

U 2 Is defined for the second sublamlnate and second ply at X 2 = 2, z = -2, 

Q 2 Is defined for the fourth sublamlnate and seventh ply at X 2 = -1 . 

PS In the first card carries the Inhomogeneous part of the condition. In sub- 
sequent cards, PS's carry the coefficients of the second, third, and fourth 
variables, respectively. NBCREGL Is an arbitrary number In the first three 
cards . 
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Data 7: NDREG, NDPLY, NDV, NDV1 , YND, ZND, STPND, NOP 


Format(4I4,3dl2.6,I4) . . .A value "2" for NDVT gives distribution of 
NDV variable of the NDPLY ply In the NDREG sublaminate with respect to 
X2-coord1nate starting at X 2 = YND computed NDP times at In Interval 
of STPND. A value "3" gives the distribution with respect to z coordinate. 
Repeat Data 7 as many times as Is required. A blank card terminates the compu- 
tations of all distributions. 

The source code Is given In appendix I. Typical Input and output for a 
ENF specimen are shown In the next two appendices. The output Is self- 
explanatory. The number of plys, real and complex roots, eigen modes and 
boundary conditions for each region are printed. The Input boundary/continuity 
conditions are also printed. 

[CBB] Is the plane stress stiffness matrix. CB11, CB12,...are Its ele- 
ments. |CB| is the vector referred to In equation 1-11 of appendix II. SNC1 
and SNC2 denote K n i and K n 2 of equation 1-12. SHM2 and SHQ2 denote 
the coefficients of Q 2 and m 2 In equation 1-9. KNI, KN2, KM1 , KM2, KQ1 , 

KQ2, KP1 , and KP2 are defined In equation 1-11. 

The amplitudes of the governing variables are given as rows In four 
columns under ZERO MODE, REGION = xx and EIGENVALUE and MODE FOR 
REGION = xx. The governing variables In each region are assumed In the order 
of U 2 , W, c> 2 2 (X 2 , z = -c) and o zz (X 2 , z = -c) for each layer In that 
region starting from top layer to the bottom layer. The different columns Indi- 
cate the mode shape. In the case of zero modes, the first, second, third, and 

0 1 2 

fourth columns denote, respectively, the coefficients of x^, x^, x^ 
and Xp. In the case of real roots only the first two columns are meaning- 
ful. They represent the coefficients of cosh(a*X2*2/H) and s1nh(a*X2*2/h) , 

In that order. If NRGN Is not "0," then the second column will also become 
zero and the first column represents the coefficient of exp(a*X2*2/H) . 

Symbols, a and H, denote the real root and the thickness of the sublaminate 
respectively. In the case of complex roots, these columns represent the coef- 
ficients of cosh(a*X2*2/H)*cos(b*X2*2/H) , cosh(a*X2*2/H)*s1n(b*X2*2/FI) , 
s1nh(a*X2*2/H)*cos(b*X2*2/H) and s1nh(a*X2*2/H)*s1n(b*X2*2/H) , respectively. 

The first and second columns would represent the coefficients of 
exp(a*X2*2/H)*cos(b*x 2 *2/H) and exp(a*x 2 *2/H)*s1n(b*X2*2/H) If NRGN Is not 
Initiated to zero. a+1*b denotes a complex root. 

The distributions are output with Y, F, Z, NR, NPLY, and NV denoting 
X2-coord1nate, value of the variable NV, z-coordlnate, sublaminate, ply, and 
the variable NV, In that order. 


APPENDIX IV: 


TYPICAL 1NU1 FOR ENF SPECIMEN 


4 2 

0.0000 0.0000 

1 

3 1 

10 0 
1 
1 

12 1 

1 

1 

1.200D-01 

0.800D-01 0.400D-04 0.600D-01 

0.600D-01 

0.600D-01 

18.8000 1.1930 1.1930 0.2600 0.2600 0.4230 0.3620 0.6230 0.6230 90.0000 

0.5000 0.5000 0.5000 0.3500 0.3500 0.3500 0.1852 0.1852 0.1852 90.0000 

0NK10 
OREAL 
3.75 
1 STORE 

1.05 
1ST0RE 

0 EXIT 

0 EIGEN 
OSTORE 
0NK10 
OCOHPLEX 
808.000000764. 

1ST0RE 

EXIT 

EIGEN 

STORE 

OPOLY 

0NK10 

OREAL 

2.5 
1ST0RE 

.59 

1ST0RE 

0 EXIT 

OEIGEN 

ODO NOT STORE 
0NK10 
OCOHPLEX 
23.200000023.3 
1 STORE 

EXIT 

OEIGEN 

ODO NOT STRORE 

OPOLY 

0NK10 

1REAL 

0NK10 

1 COMPLEX 
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OPOLY 


1 

1 

2 -. 100000 D 01 

0.000000 0.000000 

0 

1 

1 

3 -. 100000 D 01 

0.000000 0.000000 

0 

1 

1 

6 -. 100000 D 01 

0.000000 0.000000 

0 

1 

1 

10 . 100000 D 01 

0.000000 0.000000 

2 

2 

1 

1 -. 100000 D 01 

0.000000 - 1.000000 

0 

1 

1 

20 . 100000 D 01 

0.000000 0.000000 

2 

2 

1 

2 -. 100000 D 01 - 

• 300000 D -01 - 1.000000 

0 

1 

1 

30 . 100000 D 01 

0.000000 0.000000 

2 

2 

1 

3 -. 100000 D 01 

0.000000 - 1.000000 

0 

1 

1 

40 . 100000 D 01 

0.000000 0.000000 

2 

2 

1 

4 -. 100000 D 01 - 

. 300000 D -01 - 1.000000 

0 

1 

1 

S 0 . 100000 D 01 

0.000000 . 120000 D 03 

2 

2 

1 

S -. 100000 D 01 

0.000000 - 1.000000 

0 

1 

1 

60 . 100000 D 01 

0.000000 0.000000 

2 

2 

1 

6 -. 100000 D 01 - 

. 300000 D -01 - 1.000000 

0 

2 

1 

1 0.000000 

0.000000 0.000000 

3 

3 

1 

1 0.000000 

0.000000 0.000000 

0 

2 

1 

2 0.000000 

0.000000 0.000000 

3 

3 

1 

2 0.000000 

0.000000 0.000000 

0 

2 

1 

3 0.000000 

0.000000 0.000000 

3 

3 

1 

3 0.000000 

0.000000 0.000000 
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APPENDIX V 


TYPICAL OUTPUT FOR ENF SPECIMEN 
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